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Executive  Summary 


The  project  consisted  of  an  initial  three-year  effort  by  Stanford,  ACLARA,  CFD  Research  Corp., 
University  of  Utah,  and  the  University  of  California  at  Santa  Barbara;  and  an  additional  task  by 
Stanford,  Sandia  National  Labs,  and  The  Johns  Hopkins  University.  The  overall  effort  was  led 
by  Stanford  University,  funded  by  Dr.  Anantha  Krishnan  of  the  Defense  Advanced  Research 
Projects  Agency  (DARPA),  and  monitored  by  Clare  Thiem  of  the  Air  Force  Research  Laboratory 
(AFRL). 

The  goal  of  the  first  three-years  of  the  project  was  the  development  of  novel  on-chip  assay 
devices  and  modeling  capabilities  which  will  enable  optimized  design  processes  and  create  new 
methods  to  realize  robust,  field-portable  microfluidic  devices  for  the  detection.  First,  the  team 
developed,  validated  and  commercialized  new  multiphysics  models  to  the  Bio-Micro-electro- 
mechanical  systems  (MEMS)  community  (over  50  organizations)  through  CFD-ACE-I-.  These 
included  codes  on  dielectrophoresis,  electrothermal  flow,  microbead  assays,  and  surface  plasmon 
detection  devices  for  antigen/antibody  assays.  Second,  the  team  developed  rapid,  automated 
biokinetics  data  extraction  methods  for  antibody  assays  (including  lOOOx  reduction  in  time  in  a 
Biacore  Corp.  unit  as  an  example,  from  days  to  minutes).  Third,  the  team  discovered, 
formulated,  and  created  models  of  a  new  electrokinetic  instability.  This  instability  was  leveraged 
in  creating  a  rapid  micromixer  with  lOOOx  mixing  rate  over  diffusion  in  a  typical  simple  binding 
assay.  Fourth,  the  team  experimentally  quantified  and  developed  slip  velocity  models  for 
hydrophobic  channels.  Fifth,  the  team  developed  a  novel  on-chip  assay  device  that  combines 
isoelectric  focusing  and  zone  electrophoresis  to  achieve  a  2D  assay  in  l/30th  of  the  time  of  a 
traditional  system  (and  with  10,000x  less  sample).  Lastly,  the  team  developed  methods  for  on- 
chip  sample  stacking  to  improve  sensitivity.  This  effort  produced  a  fast  field  amplified  sample 
stacking  method  achieving  1  lOOx  increase  in  concentration  in  less  than  10  sec.  The  latter  effort 
was  so  successful  that  Dr.  Anantha  Krishnan  of  DARPA  provided  funding  to  extend  the  project 
for  an  additional  task  to  explore  sample  stacking  further. 

The  additional  task  focused  on  developing  rapid  sample  pre-concentration  methods  to  improve 
the  sensitivity  of  on-chip  assays  (with  a  target  figure  of  merit  of  10,000  fold  increase  by  project 
end).  The  team  developed  high  fidelity,  computational  codes  for  full  electrokinetic,  convective- 
diffusion  assays  with  multiple  species  and  fast  reaction  kinetics  capability.  Experimentally,  the 
team  developed  stacking  methods  using  field  amplified  sample  stacking,  thermal  gradient 
focusing,  and  isotachophoresis.  By  project  end,  the  team  experimentally  demonstrated  million¬ 
fold  sample  concentration  increase  (approximately  three  orders  of  magnitude  better  than  the 
second  best  ever  demonstrated  on  or  off  chip).  This  was  accomplished  using  optimized  (and 
novel  control  methods  for)  isotachophoresis.  The  team  demonstrated  and  were  able  to  inject, 
stack,  separate,  and  detect  analytes  with  100  attomolar  initial  concentrations  in  less  than  two 
minutes,  the  highest  sensitivity  ever  reported. 
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The  activities  of  the  project  can  be  categorized  into  the  areas  of  Multiphysics  Models  for 
Microfluidic  Systems,  Novel  Electrokinetic  Devices  and  Models,  and  Sample  Preconcentration 
Devices  and  Methods  for  High  Sensitivity  Assays.  The  accomplishments  are  summarized  here 
as  follows: 

1.  Multiphysics  Models  for  Microfluidic  Systems:  Developments  being  used  by  over  50 
organizations  include: 

a.  Experimentally  validated  and  commercialized  tools  for  coupled  convective- 
diffusion,  thermal,  concentration  and  electric  fields. 

b.  Models  and  theory  for  dielectrophoresis  and  electrothermal  flows. 

c.  Models  and  theory  for  sample  stacking  via  field  amplified  sample  stacking  and 
isotachophoresis. 

d.  High-fidelity  models  for  Biacore-type  biosensors  and  bimolecular  reactions. 

These  include  Biacore  sensors,  particle-based  assays,  biokinetics  data  extraction, 
and  design  for  faster,  more  sensitive  assays. 

2.  Novel  Electrokinetic  Devices  and  Models  Research: 

a.  Discovered  and  identified  a  new  electrohydrodynamic  instability  that  occurs  in 
the  electrokinetic  flow  regime  and  can  couple  with  electroosmotic  flow. 

b.  Proposed  mechanism  and  model  for  this  electrokinetic  instability  (EKI)  including 
linear  stability  analyses  and  non-linear  simulation.  Models  can  predict  critical 
conditions,  mixing  rate. 

c.  Demonstrated  robust  micromixing  devices  based  on  EKI. 

d.  Demonstrated  a  novel  miniaturized  2D  assay  combining  (for  the  first  time  on- 
chip)  isoelectric  focusing  and  capillary  electrophoresis.  Assay  uses  10  ng  sample, 
has  800  peak  capacity,  and  can  be  performed  in  ~60  min. 

3.  Sample  Preconcentration  Devices  and  Methods  for  High  Sensitivity  Assays  ResearcH: 

a.  Developed  thermal  gradient  focusing  techniques,  demonstrating  10,000-fold 
stacking  and  showing  dependence  of  stacking  on  temperature  gradient,  electric 
field,  and  other  critical  parameters. 

b.  Developed  On-chip  capillary  electrophoresis  (CE)  with  field  amplified  sample 
stacking  (PASS).  Demonstrated  12,000-fold  concentration  increase  (171x  greater 
than  any  other  on-chip  PASS).  Can  stack  and  separate  in  <20  s. 

c.  As  an  additional  task,  developed  optimized  isotachophoresis  assays.  — 
Demonstrated  sample  preconcentration  of  over  million-fold  in  less  than  2  min. 
Demonstrated  injection,  stacking,  and  detection  of  100  attomolar 
concentrations — to  the  knowledge  of  the  team,  the  highest  ever  reported 
sensitivity  of  any  electrophoretic  or  chromatographic  assay. 

The  project  had  a  very  significant  impact  on  the  microfluidics,  lab-on-a-chip  and  biosensor 
communities  as  a  whole.  In  addition  to  these  technical  achievements,  the  project  resulted  in  over 
168  journal  and  conference  papers,  three  inventions  and  patents  and  patent  applications,  and 
trained  10  PhD’s,  three  MS  graduates,  and  six  postdocs.  The  work  led  to  six  awards  from 
conferences,  national  organizations,  and  the  US  government.  In  addition,  this  project  funding 
and,  more  importantly,  the  associated  technical  accomplishments  from  this  project  contributed  to 
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the  successful  promotion  to  Associate  Professor  with  tenure  for  two  of  the  co-PIs  involved  in  the 
effort  (Myszka  and  Santiago).  Three  of  these  PhDs  and  postdocs  are  now  faculty  members  at 
major  universities  and  building  a  career  on  microfluidics  research.  Nine  others  are  engineers, 
researchers,  and/or  managers  concentrating  on  microfluidics  research  and  development  at  either 
companies  or  national  laboratories. 


1.0  Introduction 

This  final  technical  report  describes  the  major  achievements  of  the  multi-year,  multi-center, 
research  contract  entitled,  “Quantitative  Development  of  Biomolecular  Databases,  Measurement 
Methodology,  and  Comprehensive  Transport  Models  for  Bioanalytical  Microfluidics,”  funded  by 
Defense  Advanced  Research  Projects  Agency’s  (DARPA’s)  Simulation  of  Biological  Systems 
(SIMBIOSYS)  Program  (Contract  No.  F30602-00-2-0609).  The  goal  of  the  first  project  was  the 
development  of  novel  on-chip  assay  devices  and  modeling  capabilities  to  enable  optimized 
design  processes  and  create  new  methods  to  realize  robust,  field-portable  microfluidic  devices. 
Naturally  this  work  involved  a  combination  of  modeling  and  simulation,  experimentation,  and 
prototyping.  Due  to  success  of  research  in  the  area  of  on-chip  electrophoretic  sample 
preconcentration  an  additional  task  was  added  to  pursue  the  work  further.  The  body  of  the  report 
is  basically  divided  into  three  sections.  The  first  section  corresponds  roughly  to  the  original 
project,  while  the  second  and  third  sections  describe  the  experimental  and  modeling  elements  of 
the  extension. 


2.0  Quantitative  Development  of  Biomolecular  Databases,  Measurement 
Methodology,  and  Comprehensive  Transport  Models  for  Bioanalytical 
Microfluidics 

The  goal  of  the  original  project  was  the  development  of  novel  on-chip  assay  devices  and 
modeling  capabilities  which  could  enable  optimized  design  processes  and  create  new  methods  to 
realize  robust,  field-portable  microfluidic  devices  for  a  variety  of  applications.  The  project 
incorporated  an  interactive  approach  that  leveraged  experiments,  simulation,  and  fabrication  so 
that  each  step  of  the  experiment  data  acquisition,  flow  experiment  interrogations,  and  ligand- 
receptor/flow  characterization  were  modeled  in  detail  as  the  data  was  obtained.  The  goal  of  the 
project  was  successfully  accomplished  over  the  course  of  the  project,  with  devices  being 
developed  across  a  wide  range  of  applications  and  modeling  tools  developed  for  CFD-ACE+  to 
both  model  device  peformance  and  provide  guidance  for  future  design  and  development.  This 
section  briefly  reviews  some  of  the  project  team’s  accomplishements  in  this  area. 
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2.1  Development  of  Multiphysics  Models  (CFD  Research  Corporation,  U.  of  Utah,  UCSB, 
Stanford) 

The  role  of  CFD  Research  Corporation  (CFDRC)  was  to  provide  tools  to  model  the  new  devices 
and  techniques  developed  by  the  experimental  investigators.  CFDRC’s  involvement  in  the  effort 
would  also  provide  a  means  to  transition  new  tools  and  newly  found  physics  resulting  from  this 
project  to  a  broad  range  of  researchers  and  device  designers  when  incorporated  into  their  CFD- 
ACE+  modeling  software  tools.  CFDRC  developed  models  to  simulate  dielectrophoresis  (DEP) 
and  electrothermal  flow.  Both  of  these  models  have  been  used  extensively  in  the  design  and 
development  of  tunable  laser  cavity  (TEC)  sensor  (in  collaboration  with  UCSB),  design  of 
sample  stacking  system  (Stanford)  and  design  of  microfluidic  devices  by  ACLARA.  The 
summary  of  the  models  developed  and  their  validation  are  reported  below. 


2.1.1  Simulation  of  Electrothermal  flow  (CFDRC  and  UCSB) 

Use  of  microfluidic  devices  in  biological  sample  preparation  and  assays  has  received  increasing 
attention  in  recent  years.  These  devices  often  rely  on  the  fast  and  thorough  mixing  of  different 
analytes  in  a  buffer  medium.  Normally,  diffusion  is  the  mechanism  through  which  uniform 
concentrations  are  achieved.  Although  the  typical  characteristic  length  scale  associated  with 
microfluidic  devices  is  small,  diffusion  alone  cannot  provide  adequate  mixing  because  of 
inappropriate  device  dimensions,  operational  times,  and  complex  mixing  protocols.  To  reduce 
the  diffusion  length  scale,  external  forcing  is  required  to  induce  repeated  stretching  and  folding 
of  the  material  interface[2]  leading  to  both  short-term  and  long-term  chaotic  behavior.  In  large- 
scale  fluid  dynamic  systems,  turbulence  plays  a  critical  role  in  facilitating  mixing.  In 
microfluidic  devices,  however,  small  feature  size  and  lack  of  nonlinearity  make  use  of  turbulence 
impossible.  To  introduce  nonlinearities  in  the  flow,  passive  and  active  means  of  introducing 
forcing  functions  have  been  reported  in  the  literature.  For  example,  rapid  changes  in  flow 
directions  through  use  of  tortuous  flow  paths[3]  or  grooved  surfaces  [4]  has  been  tried  to  reduce 
diffusion  distances.  Stroock,  et.  al.[5],  Ajdari[6],  and  Johnson,  et.  al.[7],  have  demonstrated  a 
circulating  cellular  flow  generated  via  inhomogeneous  surface  charge  as  a  means  to  facilitate  the 
mixing  process  using  Electroosmosis.  Oddy,  et.  al.[8]  used  conductivity  gradients  in  the  bulk 
fluid  as  a  method  to  induce  chaos  and  mixing  when  an  alternating  electric  field  is  applied. 
Magnetic  fields  [9]  have  also  been  used  to  facilitate  mixing.  Other  active  means  include  use  of 
bubbles  (thermally  actuated  [10]  and  more  recently  gaseous  slugs  [11])  and  acoustic  stirring[12]. 
Recently  it  has  been  shown  that  use  of  a  motile  bacterium  can  enhance  mixing  in  microfluidic 
flow[13].  A  comprehensive  review  of  mixing  in  microfluidic  devices  can  be  found  in  reference 
[14]. 


When  an  AC  electric  field  is  applied  to  an  electrolyte  solution,  several  modes  of  AC 
electrokinetic  phenomena  can  occur,  all  of  which  primarily  depend  on  the  applied  voltage  and 
frequency.  When  the  applied  voltage  is  low  and  the  frequency  is  sufficiently  high, 
electrothermal  flow  dominates  the  electrolysis  and  AC  electroosmosis[15].  Ramos  et.  al.[16] 
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formulated  and  discussed  in  detail  the  electrothermal  flow  induced  in  the  proximity  of  an 
electrode  that  is  under  the  influence  of  the  AC  electric  field.  When  it  is  applied  to  an  electrolyte 
solution,  the  solution  is  subjected  to  coulomb ic  and  dielectric  forces  due  to  variations  in  the 
dielectric  properties  (electrical  permittivity  and  conductivity)  as  a  function  of  temperature. 
Though  joule  effect  is  often  the  primary  mechanism  for  heating  the  solution,  other  external 
sources,  such  as  laser  illumination,  can  induce  thermal  gradients,  and  hence,  electrothermal 
flows.  Electrothermally  induced  flow  depends  on  a  variety  of  parameters  such  as  electrical 
conductivity  and  permittivity,  the  voltage  and  frequency  of  the  applied  electric  field,  and 
electrode  configuration  and  channel  dimensions. 


Mathematical  formulation 


Let’s  consider  a  microfluidic  channel  consisting  of  a  pair  of  co-planar  electrodes  or  electrode 
arrays  kept  at  the  bottom  floor  of  the  channel  as  shown  in  Figure  2.1.  The  electrodes  have  a 
width  of  d2  and  a  gap  of  2di  between  adjacent  electrodes.  They  are  energized  by  an  electric 
potential  of  -Vo  and  Vo,  respectively.  The  force  density  due  to  application  of  the  electric  field  is 
given  by [17]: 

in  1  j 

/  =  /7E — EVf  +  -V 
“2  2 


Pm 

E' 

UcJ 

T 

(2.1) 


Where  pe  is  the  charge  density,  E  is  the  electric  field,  s  is  the  permittivity  of  the  buffer,  and  pm  is 
the  buffer  density.  The  first,  second,  and  third  terms  on  the  right  side  of  the  above  equation 
represents  the  Coulomb  force,  the  dielectrophoretic  force,  and  the  electrostriction  force, 
respectively.  If  the  applied  electric  field  is  an  alternating  one,  then  the  averaged  force  density 
can  be  expressed  as[16]: 


<7e{a-  P) 
a  +  ims 


(VT-E„)E„--m|Ej  Vr 


(2.2) 


Linear  variations  of  the  electric  properties  are  assumed  in  deriving  the  above  equation.  In 
Equation  2,  a  is  electrical  conductivity,  T  is  the  temperature,  and  a  and  P  are  the  coefficients  of 
variations  of  the  electrical  properties  with  respect  to  the  temperature.  The  coefficients  a  and  p 
are  expressed  as: 


\  ds  ^1  d(7 

a  = - ,  P  = - 

sdT  (7  dT 


(2.3) 


This  report  focuses  only  on  joule  heating  (/=aE^)  as  a  source  of  introducing  thermal  gradients. 
The  conservation  equation  for  thermal  field  is  expressed  as: 


dT  2  2 

p  c - h  yC>  c  (u  •  Vjr  =  kV  T  +  crE 

I  m  p  ^  /  m  p  ^  ^ 

ot 
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Where  Cp  is  the  specific  heat  capacity,  u  is  the  velocity  vector,  and  k  is  the  thermal  conductivity. 
The  fluid  flow  is  governed  by  Navier-Stokes  and  continuity  equations: 


du 

dt 


+  (u-V)u 


=  —Vp  +  //V^u  +  (f ) 


Vu  =  0 


(2.5) 


Where  p  is  pressure  and  <f>  is  the  momentum  source  induced  due  to  the  electrothermal  flow. 
The  electric  field  satisfies  the  following  relations: 

dp 

E  =  -VF,-^  +  V-j  =  0 

St  (2.6) 

V  •  D  =  0,  D  =  sE,  j  =  (jE 


Where  v  is  the  electric  potential,  j  is  the  current  density,  D  is  the  displacement  vector,  and  a  is 
the  electrical  conductivity  of  the  buffer.  The  solution  to  the  coupled  Equations  (4)  through  (6), 
subjected  to  appropriate  boundary  conditions,  describe  the  thermal,  electric,  and  flow  fields. 
Two  modeling  approaches  are  used  to  characterize  the  mixing  of  analytes  in  microfluidic 
devices.  When  the  size  of  the  analyte  is  very  small  (sub-micron)  and  the  inertial  acceleration 
relative  to  the  fluid  can  be  neglected,  a  continuum-based  advection/diffusion  equation  for  the 
analyte  concentration  is  solved.  When  the  analyte  molecules  are  larger  (micron-sized)  and  the 
inertial  acceleration  is  important,  a  kinematic  equation  describes  the  motion  of  the  particles  [2] 


Exact  solution 

Equations  (2.4)  through  (2.6)  are  coupled,  non-linear,  PDEs  that  usually  preclude  analytical 
solution.  However,  this  description  follows  the  procedure  used  by  Ramos,  et.  al.  [16],  and 
rewrites  the  equations  for  low  Reynolds  number  («1)  steady- state  flow  commonly  seen  in 
microfluidic  devices.  The  resulting  linear  equations  admit  an  analytical  solution  for  simplified 
geometries.  In  this  section,  the  analytical  solutions  for  a  co-planar  electrode  structure  as  shown 
in  Figure  2.1  are  briefly  presented.  The  normalized  electric  potential  ^(normalized  by  Vo) 
distribution  is  given  by  the  Fourier  integral: 


<f>{x,y)  =  [  A(A)A.  ‘e  sin(/ix)dk 

J  0 


where. 


and 
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Here  and  Jn  are  the  Legendre  and  Bessel  functions  of  order  n,  K  is  the  complete  elliptic 
function  of  the  second  kind.  The  general  solution  to  temperature  is  expressed  as: 


a(f)^ 

2k 


+  Ax  +  By  +  C 


(2.10) 


where  A,  B,  and  C  are  constants  determined  from  thermal  boundary  conditions.  By  assuming  an 
isothermal  condition  on  the  electrode  surface,  and  solving  heat  conduction  equation,  it  becomes: 


T  =  Z 


Ik 


(2.11) 


The  solution  for  the  electrothermal  force  is  calculated  from  the  Equation  (2.2)  by  the  Green’s 
function  method 


^.(■^.J)  =  JJ, 


G  (x,x. 

s,  ^ 


3^0,2  )/,  K.P  3^0,2 
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(2.12) 


Where  Gij  is  green’s  function  or  Stokeslet  in  the  presence  of  a  planar  boundary.  The  above 
methodology  can  be  extend  to  a  periodic  array  of  electrodes  and  calculate  the  electric  potential  in 
terms  of  Fourier  series: 


(^(x,y)  =  X  — ^cos/l  e  =n-  — 

2  "  ”  2 


(2.13) 


Where, 

P  (cos  c) 

<=—  ^  (2.14) 

k:(cos-) 

2 

Similarly,  the  flow  field  is  calculated  by  (2.12)  in  which  the  two-dimensional  singly  periodic 
Stokeslet  in  the  presence  of  a  planar  boundary.  The  above  analytical  solutions  are  used  to 
validate  the  numerical  method. 


Numerical  method 


Microfluidic  devices  usually  have  complex  electrode  patterns,  which  renders  use  of  the 
analytical  solutions  derived  based  on  the  linearization  assumption  invalid.  In  this  scenario,  the 
full-set  of  governing  equations  needs  to  be  solved  using  numerical  methods.  A  finite-volume 
method-based  solver  available  in  the  commercial  software  CFD-ACE+[18]  is  used  to  solve  the 
Equations  (2.4)  through  (2.6)  in  a  coupled  manner.  Analyte  mixing  is  studied  using  an 
advection-diffusion  equation  or  kinematic  equations.  The  basic  flow,  electric,  and  thermal 
solvers  are  available  in  CFD-ACE-l-.  The  source  term  to  the  momentum  equations  corresponding 
to  the  electrothermal  flow  has  been  added.  The  numerical  method  has  been  validated  with  the 
help  of  the  analytical  solution  described  in  the  previous  section. 
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Model  validation 


A  schematic  of  the  co-planar  electrode  array  design  is  shown  in  Figure  2.1.  Analytical  solutions 
of  this  configuration  were  used  to  validate  the  full-scale  numerical  simulation.  The  electrode 
width  and  the  gap  used  in  the  study  are  50pm.  A  DC  signal  of  ±2V  is  applied  to  the  adjacent 
electrodes.  The  ambient  temperature  is  assumed  to  be  300K.  In  all  the  simulations,  the  KCl 
solution  is  used  as  the  buffer.  It  has  a  density  of  1000  kg/m^,  viscosity  of  10'^  m^/s,  thermal 
conductivity  of  0.6  W/m-K,  heat  capacity  of  4100  J/kg-K,  electrical  conductivity  of  0.056/0hm- 
m,  relative  permittivity  of  78,  a  value  of  -0.004  /K  and  |3  value  of  0.02  /K.  Comparison  between 
numerical  and  analytical  solution  in  terms  of  velocity  and  temperature  along  a  plane  at  an 
elevation  of  5  pm  and  25  pm  from  electrode  surface  is  shown  in  Figure  2.2.  In  numerical 
simulation,  the  grid  used  is  400  by  200,  and  further  increase  of  grid  numbers  does  not  improve 
the  results.  The  numerical  and  analytical  results  agree  with  the  error  less  than  0.1%.  A  pair  of 
circulating  vortices  is  observed  very  close  to  the  electrode  edges.  The  maximum  velocity  in  this 
region  is  approximately  4  pm/s.  The  maximum  temperature  difference  is  very  small  (0.2  K). 


Figure  2.1:  Coplanar  electrode  used  in  the 
validation  study. 

Solid  lines  illustrate  electric  field;  dashed 
lines  illustrate  streamlines. 
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Figure  2.2:  Comparison  of  numerical  results  (lines)  with  analytical  data  (circles). 

(a)  Temperature  field  over  eleetrode  design  of  Figure  2.1.  (b)  Velocity  field  in  the  same  region. 


2.1.2  Model  for  Dielectrophoresis 

The  team  has  accomplished  following  in  this  task: 

•  Models  to  simulate  conventional  and  traveling  wave  DEP  have  been  implemented.  For 
this  purpose,  effective  complex  permittivity  has  been  used. 

•  The  “electrical  module”  in  CFD-ACE+  has  been  enhanced  to  predict  high  frequency  AC 
electrical  fields  accurately. 

•  The  module  was  first  tested  and  validated  by  comparing  numerical  solution  with 
analytical  results  for  a  simple  DEP  configuration.  Following  this,  simulations  of  UCSB 
Tunable  Laser  Sensor  was  performed  for  design  optimization  studies. 


DEP  phenomenon,  in  contrast  to  electrophoresis  of  charged  particles  in  electric  field,  describes 
the  motion  of  polarizable  particles  under  the  influence  of  an  non-uniform  electric  field  [19].  It 
varies  significantly  with  applied  electric  potential  and  frequency,  as  well  as  dielectrical 
properties  of  the  particles/cells  and  the  medium.  Under  the  action  of  an  externally  imposed  AC 
electric  field,  the  in-phase  component  of  the  dipole  moment  induces  conventional  DEP  (c-DEP) 
forces  while  the  out-of-phase  component  induces  traveling-wave  DEP  (tw-DEP)  forces  [20,  21]. 
The  c-DEP  and  tw-DEP  microchips  show  great  potential  for  characterizing  and  sorting  cells  and 
microorganisms  [22-24].  Standard  microfabrication  techniques  and  materials  are  used  in 
fabrication  of  DEP  devices.  Compared  with  electrophoretic  (or  electroosmotic)  separation 
systems,  DEP  devices  operate  using  low  voltage  AC  rather  than  high  direct  current  (DC)  voltage. 
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Furthermore,  DEP  devices  can  be  integrated  with  other  microfluidic  components  to  provide  an 
efficient,  rapid,  and  automatic  analysis  system. 

Design  of  DEP  devices  poses  several  challenges  because  DEP  forces  are  a  complex  function  of 
applied  electric  field  and  dielectric  properties  of  the  buffer  and  particles.  The  study  of  electric 
fields  generated  by  the  electrode  structure  becomes  essential,  particularly  in  guiding  the  design 
of  microelectrode  array.  Existing  analytical  and  numerical  models  focused  on  establishing  the 
significance  of  microelectrode  on  the  equilibrium  position  and  transient  motion  of  particle 
mixture  under  various  experimental  conditions.  They  illustrate  spatial  variations  of  DEP  forces 
within  a  device,  but  are  limited  to  specific  geometries  or  approximate  boundary  conditions  [25, 
26].  As  the  complexity  of  the  devices  increases,  practical  applications  of  these  models  are 
limited  to  predicting  levitation  height,  particle  velocity,  and  selected  particle  trajectories  [27]. 
On  the  other  hand,  high  fidelity  simulations  provide  deeper  insight  into  various  competing  forces 
and  facilitate  the  selection  of  optimal  geometry  and  process  parameters  for  maximum 
performance. 

For  thorough  elaboration  of  DEP  theory,  see  Pohl  [19];  this  section  offers  a  necessarily  brief 
account  of  theoretical  fundamentals  relevant  in  the  context  of  this  work.  The  particles  or  cells 
experience  DEP  force  from  the  non-uniform  electric  fleld  arising  due  to  geometry  or  electric 
field  non-uniformity.  The  DEP  force  can  be  derived  by  either  the  effective  dipole  moment 
method  or  the  Maxwell  stress  tensor  method.  In  general,  the  time-averaged  DEP  force  for  a 
spherical  particle  in  AC  electric  field  is  given  as  follows  [28]: 

f=L(m-V)E  (2.15) 

where  m  is  the  dipole  moment  and  E  is  the  electric  field.  For  an  isotropic,  homogeneous 
dielectric  spherical  particle,  the  DEP  force  in  general  is  expressed  as  [29]: 

fo,p  =  27ra^s„So  [o.5Re(^)V(E^o  +  +  E^o)  +  Im(ir)(EjoV(p,  +  E^^^Vcp^  +  E^^^Vcp,)]  (2.16) 

where  a  is  the  radius  of  the  particle,  8  is  the  permittivity;  subscript  m  represents  the  media,  sq  is 
the  permittivity  of  free  space,  (E^Q,EyQ,E^Q)  and  ((p^,(Py,(p^)are  the  magnitudes  and  phases  of  x, 

y,  and  z  components  of  the  electric  field;  V  is  the  gradient  operator;  and  K  is  Clausius-Mossotti 
factor.  Equation  (2.16)  shows  two  terms  to  the  DEP  force  vector:  a  real  part  related  to  the  in- 
phase  component  of  the  induced  dipole  moment  and  an  imaginary  part  related  to  the  out-of-phase 
component.  The  in-phase  component  directs  the  particles  to  the  regions  of  electric  field  maxima 
(Re(iO>0)  or  minima  (Re(iQ<0).  The  motion  of  the  particles  due  to  this  component  is  termed 
conventional  DEP  (c-DEP).  The  out-of-phase  component  directs  the  particles  to  the  regions  of 
where  the  phases  of  the  field  component  are  larger  (Im(.^  >  0)  or  smaller  (lm(^<0).  The  linear 
motion  of  the  particles  due  to  this  component  is  termed  traveling- wave  DEP  (tw-DEP).  The 
Clausius-Mossotti  factor  (K)  is  expressed  as: 
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(2.17) 


where  s*=  8  -  /cr/co,  is  the  complex  permittivity;  i  =  4^i ,  a  is  the  conductivity,  and  co  is  the 
angular  frequency  of  the  applied  field.  A  typical  Clausius-Mossotti  spectrum  for  a  microbead 
suspended  in  an  aqueous  buffer  with  Sp  =  40,  8m  =  80,  ap  =  l(S/m)  and  CTm  =  0.001  (S/m)  is 
shown  in  Figure  2.3.  At  low  frequencies,  these  beads  are  generally  more  polarizable  than  the 
suspending  medium  and  experience  positive  c-DEP.  As  the  electric  field  frequency  increases, 
the  c-DEP  force  acting  on  the  particles  changes  from  positive  to  negative.  At  the  crossover 
frequency,  when  the  real  part  of  Clausius-Mossotti  factor  is  zero,  beads  experience  no  DEP 
force.  This  frequency  critically  depends  on  the  dielectric  properties  of  the  particle  and  the 
medium.  The  c-DEP  force  has  two  components:  a  normal  component  that  levitates  (attracts)  the 
particle  in  a  direction  normal  to  the  electrode  surfaces  and  a  horizontal  component  that  pushes 
(pulls)  the  particle  away  from  (toward)  electrodes.  Both  components  of  c-DEP  forces  decrease 
significantly  as  particles  move  away  from  the  electrode.  On  the  other  hand,  the  tw-DEP 
spectrum  exhibits  a  unique  behavior.  The  imaginary  part  vanishes  at  small  and  high  values  of 
frequencies.  For  a  given  particle-buffer  pair,  the  imaginary  part  is  maximum  only  at  one  critical 
frequency.  In  the  case  of  biological  cells,  one  may  assume  cell  as  an  ohmic  spherical  particle 
and  neglect  the  trans-membrane  conductance  and  surface  conductivity.  This  assumption  is 
particularly  good  for  mammalian  cells  that  do  not  have  a  cell  wall  and  have  a  membrane 
thickness  at  least  three  orders  of  magnitude  smaller  than  the  cell  radius. 


-0.40 


0.80 


0.40 


0.00 


Frequency  (Hz) 


Figure  2.3:  Clausius-Mossotti  spectrum  for  a 
microparticle  with  8p  =  40,  Sm  =  80,  ap  =  l(S/m) 
and  Qm  =  0.001  (S/m). 


Figure  2.4:  Illustration  of  various  forces  acting  on  a 
particle  kept  near  an  array  of  electrodes  in  a  channel. 
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Comparison  of  Force  Calculation  with  Analytical  Solution  in  an  Interdigitated  Electrode 
Structure 


Figure  2.4  illustrates  various  forces  acting  on  a  particle  in  a  microchannel  device.  The  channel 
has  an  array  of  interdigitated  electrodes  at  the  bottom.  An  AC  electric  field  is  applied  with  1 80 
degrees  shifted  in  phase  between  consecutive  electrodes.  The  polarizable  particles  levitate  above 
the  electrode  under  a  negative  DEP  force  Fy.  This  force  is  balanced  by  a  gravitational  force  Fg^ 
which,  for  a  spherical  particle  of  radius  a,  is  Fg  =  3/4  %a  Apg,  where  Ap  is  the  difference  between 
the  density  of  the  particle  and  the  suspending  solution,  and  g  is  the  gravitational  acceleration.  In 
a  traveling  electric  field,  the  levitated  particles  can  also  move  horizontally  along  the  electrode 
array.  Particles  of  different  polarizability  move  along  the  electrode  with  different  velocities. 
Under  steady-state  conditions,  this  force  is  balanced  by  Stokes  drag:  Fstoke  =  6Tir\av,  so  that 
Fx+Fstoke  =  0  where  v  is  the  traveling  wave  velocity  and  p  is  the  viscosity  of  the  suspending 
medium.  Simulation  of  dynamic  motion  of  a  particle  involves  balancing  these  forces  and 
computing  the  location  of  the  particle.  The  channel  height  is  sufficiently  large  when  compared 
to  electrode  width  and  spacing,  so  that  the  boundary  condition  at  infinity  is  used  at  the  top 
channel  wall.  The  peak-to-peak  potential  is  2  Vq.  An  assumption  is  that  the  thermal  conductivity 
of  the  metal  electrode  is  much  larger  than  that  of  the  fluid  so  that  the  temperature  deviates  very 
little  from  the  electrode  temperature  Ts,  The  temperature  deviation  in  the  fluid  due  to  heating  of 
electrical  current  is  of  interest.  Both  electrical  and  thermal  field  are  analytically  solved  and  the 
detailed  are  published  elsewhere  [27]. 
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Figure  2.5:  Comparison  of  analytical  and  numerical  results  for  the  normalized  potential,  x-component 

electrical  field  and  DEP  forces  on  a  sphere. 


Figure  2.5  shows  comparison  of  analytical  and  CFD-ACE+  (numerical)  results  and  they  agree 
very  well.  It  is  also  observed  that  the  CFD-ACE+  is  able  to  capture  the  singular  behavior  of  the 
electrical  field  and  the  DEP  force  near  the  electrode  edge,  where  they  usually  tend  to  go  to 
infinity.  The  simulation  also  predicted  the  levitation  force  fy  to  be  symmetric  while  the  lateral 
force  fx  to  be  anti-symmetric,  with  respect  to  the  edge  and  the  center  of  the  electrode  when  the 
electrode  width  is  equal  to  that  of  the  spacing.  This  behavior  has  not  been  predicted  by 
approximate  solutions  reported  earlier  by  other  investigators. 


2.1.3  Model  and  Investigation  of  Slip  Condition  on  Hydrophobic  Surface 

Nano-bubbles  have  recently  been  observed  experimentally  on  smooth  hydrophobic  surfaces, 
cracks  on  a  surface  can  likewise  be  the  site  of  bubbles  when  partially  wetting  fluids  are  used. 
Because  these  bubbles  may  provide  reduced  stress  boundary  condition  and  modify  considerably 
the  friction  generated  by  the  solid  boundary,  it  is  of  interest  to  quantify  their  influence  on  fluid 
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flows  over  these  slip  boundaries.  This  nano-bubble  theory  has  reeently  been  adopted  to  interpret 
experiments  of  flow  over  a  hydrophobie  channel  [30], 


A  schematic  of  the  numerical  treatment  of  the  slip  problem  is  shown  in  Figure  2.6.  Because  of 
the  flow  surface  interaction,  Stokes  slip  boundary  condition  is  used  to  replace  the  traditional  no¬ 
slip  condition  on  the  surface  modified  walls,  particularly  for  flow  in  micro-sized  channels  [31]. 


(2.18) 


Here  the  coefficient  X  is  termed  as  slip  length.  Normal  to  the  surface,  the  non-penetration 
condition  still  applies,  i.e. 

=  0  (2.19) 


1 


Boundary  conditions  can  be  converted  into  a  discrete  form  at  the  boundary  cell  center  by  finite 
difference.  Writing  n  as  the  normal  vector,  ti  and  xi  as  two  principal  tangential  vectors  on  the 
surface,  the  two  components  of  slip  condition  becomes: 


1  1 
r  ,  u„  -u 


u^-r,=A- 


8 

r  r 

U^-T2=  a  _  -Tr, 


5 


(2.20) 


where  the  subscript  c  represents  for  cell  center  and  w  for  wall.  These  equations  along  with 


1 


n=0 


(2.21) 


determine  the  wall  velocity  components  in  terms  of  those  at  the  cell  center. 
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Figure  2.6:  Schematic  of  numerical  treatment  of  the  slip  condition  on  surface. 


For  2D  flow,  the  relation  is: 


u  = 


-Xn^[-nyU^+n^v^) 


X  +  6 


V.,  = 


X  +  d 


(2.22) 


Validation  Case: 


A  simple  channel  flow  is  considered,  in  which  slip  occurs  on  top  of  the  channel  wall.  This 
problem  admits  analytic  solution: 


_  \  dp  h  +  2X 

2p  dx\_  h  +  X 


hy-y^ 


(2.23) 


Simulation  results  of  dimensionless  velocity  distribution  agree  excellently  with  analytic  solution. 
The  slip  velocity  on  the  top  wall  is  given  by: 


dp  X 
2p  dx\_h  +  X 


For  small  slip  length,  the  linear  and  quadratic  approximation  gives: 


Xh  dp 
2p  dx  ’ 


Xh  dp  X 
2p  dx\_  h 


(2.24) 


(2.25) 
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These  results  are  plotted  in  Figure  2.7. 


-(|j/h^)/(dp/dx)  Slip  Length  X 

Figure  2.7:  Comparison  of  normalized  velocity  across  the  channel  and  the  slip  velocity  on  top  channel 

wall. 


2.1.4  Simulation  of  Tunable  Laser  Cavity  Sensor  (CFDRC  and  UCSB) 


Problem  Description 

Immunoassays,  which  rely  on  specific  antigen-antibody  binding  for  identification  of  proteins  in  a 
sample,  have  applications  in  both  clinical  laboratories  for  medical  diagnostics  and  treatment 
monitoring,  and  in  research  laboratories  for  highly  multiplexed  testing,  such  as  for  biomarker 
identification.  In  these  cases,  throughput  is  a  key  consideration.  One  factor  limiting  throughput 
is  diffusion  of  analyte  to  the  reporter.  An  incubation  step  of  minutes  to  hours  allows  diffusion- 
limited  binding  to  reach  detectable  levels.  These  tests  are  usually  performed  at  centralized  labs 
where  high  throughput  is  achieved  through  robotics  and  highly  parallel  assays.  However,  if  the 
assay  could  be  moved  from  the  centralized  lab  to  the  point  of  care  -  realizing  immunoassay- 
based  diagnostics  “on  demand”  -  the  test  must  be  faster,  as  well  as  smaller,  while  maintaining 
high  sensitivity. 

In  response  to  this  need,  Microfluidic  assays  for  diagnostics  have  developed  dramatically  in 
recent  years,  allowing  the  realization  of  the  lab-on-a-chip  concepts  for  point-of-care  diagnosis 
and  high  throughput  screening  for  molecular  diagnostics.  The  tunable  laser  cavity  (TLC)  sensor, 
being  developed  at  UC  Santa  Barbara,  uses  a  semiconductor  laser  coated  with  a  layer  of  capture 
antibodies  to  detect  the  presence  of  antigen  in  the  flow.  When  the  antigen  antibody  complex  is 
formed,  the  laser  frequency  changes  with  respect  to  the  reference  frequency.  In  order  to  enhance 
the  sensitivity  of  the  sensor,  it  is  desired  to  increase  the  concentration  of  antigen  near  the  binding 
site  where  antibody  is  coated.  In  this  regard,  CFDRC  modeled  UCSB  system  to  understand  the 
effect  of  DEP  and  electrothermal  forces  on  the  assay  performance. 
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This  report  presents  the  simulation  results  using  DEP  model.  The  simulation  parameters  are 
listed  below.  In  the  first  set  of  simulation,  velocity  field  and  particle  trajectories  are  computed 
for  an  idealized  system  that  consists  of  an  expansion-contraction  nozzle.  A  peak-to-peak  voltage 
of  20  volts  at  the  frequency  of  lOkHz  is  applied.  Other  parameters  are  listed  in  Table  2.1.  The 
team  randomly  injects  hundred  particles  close  to  the  throat  of  the  channel  and  solve  for  their 
location  as  time  evolves.  The  induced  fluid  flow  is  calculated  by  solving  the  Navier-Stokes 
equation. 


Table  2.1:  Physical  parameters  used  in  the  UCSB  TLC  Sensor  study. 


Ps 

a 

P 

O'™ 

«0 

1000 

1.5E-6 

lE-3 

0.056 

lE-18 

71/8 

80 

2.61 

Results  for  particle  configuration  at  four  different  time  instants  are  shown  in  Figure  2.8.  The 
induced  flow  pattern  is  also  shown.  Since  the  particle  conductivity  is  smaller  than  that  of  the 
media,  the  DEP  force  alone  will  drive  particles  to  regions  of  smaller  electrical  field  (negative 
DEP),  as  indicated  in  Figure  2.8.  The  magnitude  of  fluid  flow  is  about  50  pm/s. 
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Figure  2.8;  Particle  position  at  various  instants  in  time  throughout  the  first  40  ms  of  the  DEP  process. 

The  maximum  velocity  is  approximately  50  pm/s. 

These  results  have  been  validated  by  UCSB  whose  Particle  Image  Velocimetry  (PIV) 
measurements  predicted  a  fluid  velocity  of  approximately  60  pm/S  (However,  it  was  pointed  out 
that  the  numerical  simulation  of  by  UCSB  computed  a  fluid  velocity  of  70  pm/s  if  electrothermal 
effects  were  included,  and  hence  the  velocity  due  to  particle  dielectrophoresis  was  therefore 
estimated  as  10  pm/s  [2]).  In  literature,  similar  experiment  by  Prof  Green’s  group  (2000), 
calculated  that  the  electrothermal  effects  due  to  light  illumination,  resulted  in  velocities  of 
approximately  10  pm/s.  Because  the  variation  in  the  estimated  velocities  was  significant,  it  was 
decided  that  the  team  include  models  for  electrothermal  effect  in  the  simulation  and  this  is 
explained  in  the  next  section. 


2.1.5  Analysis  and  Modeling  of  UCSB  System  for  Electrothermal  Flow  Effects 

This  section  presents  the  results  from  the  electrothermal  flow  simulations  in  Figure  2.9,  where 
the  electric  potential  and  temperature  distributions  and  the  velocity  vector  have  been  plotted. 
The  AC  frequency  of  10^  Hz  and  the  RMS  voltage  of  5V,  which  corresponds  to  a  7.07  V  peak 
voltage,  were  used  in  the  simulation. 
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Figure  2.9:  Computed  potential,  temperature,  and  velocity  field  data  for  the  (|)  =5 

Vrms  case. 

The  velocity  and  electric  field  of  the  electrothermal  flow  result  in  increased  binding 

rate. 

It  is  noted  that  the  velocity  field  near  the  electrodes  is  diverted  toward  the  binding  site  where  the 
antibody  is  coated.  As  a  result,  enhanced  surface  binding  occurs  when  compared  to  (cross¬ 
stream)  diffusion  driven  binding  in  a  static/uni-direction  flow  environment.  A  closer  view  of  the 
velocity  field  is  shown  in  Figure  2.10  to  illustrate  the  bending  of  the  streamlines.  Since  the 
channel  flow  is  superimposed  on  the  electrothermal  flow,  the  team  recommended  that  the 
binding  site  may  be  moved  the  downstream  where  flow  convection  was  more  pronounced 
(thereby  reducing  mass  transport  limitations). 
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Electrode  Binding  Site  Electrode 


Figure  2.10:  A  closer  view  of  velocity  field  near  the  binding  site  (same  conditions  as  Figure  2.9). 


The  time  evolution  of  surface  concentration  of  antigen-antibody  complex  is  plotted  in  Figure 
2.11  and  compared  to  the  case  where  no  electrothermal  flow  is  involved.  For  small  periods,  the 
increase  of  surface  concentration  is  very  significant.  The  electrothermal  flow  enhances  binding 
significantly.  For  example,  the  time  required  to  achieve  90%  of  maximum  surface  concentration 
in  the  presence  of  electrothermal  flow  is  25%  shorter  than  that  without  electrothermal  flow. 
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Figure  2.11:  Comparison  of  surface  concentration  of  antigen-antibody  for  (a)  Sshort  and  (b) 

long  time  periods. 


The  surface  concentration  of  antibody  is  1 .7E-3  nM/cm. 


Figure  2.1 1  shows  the  concentration  contour  of  the  antigen  in  the  charmel.  Note  that  the  antigen 
flows  toward  the  binding  site  because  of  the  electrothermal  effect.  This  is  in  sharp  contrast  to 
how  the  antigen  is  transported  in  the  case  when  there  is  no  electrothermal  flow.  The  local 
mixing  enhancement  helps  to  eliminate  mass  transport  limitations  leading  to  increased  surface 
coverage.  The  higher  the  electrical  potential  is,  the  more  pronounced  is  the  surface  binding. 
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From  the  numerical  simulation  studies,  the  team  concluded  that  the  electrothermal  flow  is  more 
effective  in  enhancing  surface  binding  in  a  TLC  sensor.  The  net  effect  is  reduced  detection  time. 
By  increasing  the  applied  voltage,  it  is  possible  to  accelerate  the  binding  further.  The  current 
simulations  clearly  indicate  that  the  DEP  forces  are  not  significant  in  this  sensor  due  to  the  size 
of  the  antigen  molecules. 

2.1.6  Simulation  of  Field  Amplified  Sample  Stacking 

Field  amplified  sample  stacking  (FASS)  is  an  important  technique  in  analytical  chemistry, 
especially  in  modem  microfluidic  technology.  Recent  development  in  biochip  industry  focuses 
on  manipulation  of  extremely  small  amount  of  sample  in  integrated  microfluidic  devices.  One  of 
the  challenges  concerning  this  technology  is  how  to  obtain  a  measurable  signal  from  such  very 
low  concentration  of  sample.  To  enhance  the  sensitivity  of  the  assay,  it  is  cmcial  to 
preconcentrate  the  sample  before  conducting  the  experiments.  FASS  is  one  such  technique 
widely  used  on  the  microfluidic  platform. 

The  physical  mechanism  for  the  FASS  lies  in  variations  in  the  electrical  field  in  the  different 
regions  of  the  microchannel.  The  process  is  shown  schematically  in  Figure  2.12.  The  channel 
contains  two  types  of  buffer  with  a  large  difference  in  electrical  conductivity.  The  buffer  in  the 
middle  of  the  channel  has  lower  conductivity  and  the  sample  is  injected  in  this  region,  while  the 
buffer  in  other  parts  of  the  channel  has  a  higher  electrical  conductivity.  Once  the  electrical 
potential  is  applied  across  the  channel,  the  electrical  field  inside  the  channel  varies  significantly 
due  to  variations  in  the  electrical  conductivity.  The  electrical  field  is  high  in  the  sample  region 
and  therefore,  the  sample  moves  faster  toward  the  boundary  between  the  sample  and  the  mnning 
buffer  of  higher  conductivity.  As  the  sample  accumulates  at  the  interface,  the  electrical  field 
decreases  and  the  sample  migration  slows  down.  Moreover,  the  electroosmotic  flow,  which  is 
directly  proportional  to  the  electrical  field  will  also  become  weak,  causing  further  slowing  down 
of  the  sample  migration.  As  a  result,  the  sample  will  stack  into  a  thin  band  near  the  interface 
between  the  buffer  and  sample  region.  The  exact  location  (forward  or  rear  interface)  depends  on 
the  sample  charge.  The  concentration  of  sample  in  this  thin  band  is  much  higher  compared  to 
that  before  the  electrical  field  is  applied. 
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Figure  2.12:  Illustration  of  PASS  in  a  microfluidic  channel. 


Governing  Equations  Relevant  to  Stacking  Phenomenon 


The  PASS  can  be  described  by  the  electrokinetic  theory.  The  transport  of  the  sample  and  running 
buffers  obeys  mass  conservation  laws: 


- C:  +- 

dt 


dx; 


-J.  =0 


(2.26) 


where  the  flux  (Jy)  of  species  i  is  given  by  Nemst-Planck  equation: 


0^  dc. 

J..  =u  c.  -  z.c.co - D  — - 
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OX.  ox . 


(2.27) 


Here  c,,  A,  and  <»/  are  the  concentration,  difflisivity  and  electrophoretic  mobility  of  the  specie  i, 
respectively,  u  is  the  fluid  velocity  and  (j)  is  the  electrical  potential.  In  the  above  equation,  it  is 
assumed  that  the  double  layer  at  the  particle-electrolyte  surface  is  thin,  the  electrolyte  is  fully 
ionized  and  the  solution  is  dilute.  Consequently,  the  introduction  of  charged  species  will  not 
alter  the  properties  such  as  density,  viscosity  etc. 

The  motion  of  background  electrolyte  is  described  by  the  Navier- Stokes  equation  for  the 
incompressible  Newtonian  fluid  flow: 


d  d  dp  d 

—  pu.-\ - pu  u.  = - 1 - p 

dt  '  dx.  ^  '  dx.  dx. 

J  I  I 


du.  du. 
— =-  +  — - 

dx .  dx. 


2  dip 
+  SK  g  — 

dx. 


(2.28) 
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where  p  is  the  density,  p  is  the  hydrodynamic  pressure,  k~'  is  the  Debye-Huckel  thickness,  s  is 
the  dielectrical  constant  of  the  medium  and  ^  is  the  zeta  potential.  The  last  term  on  the  right 
hand  of  equation  is  the  electroosmotic  volumetric  source  term  arises  due  to  the  interaction  of 
charged  wall  and  the  ions  in  the  solution.  Two  different  approaches  to  simulate  PASS  have  been 
used. 


Approach  1  (Electroneutrality  Assumed) 


In  the  first  approach,  electroneutrality  assumption  has  been  used  and  solved  for  the  current 
conservation  equation  or  conduction  equation.  The  governing  equations  are: 

V-z^  =  0  (2.29) 

where  i  is  the  electrical  current  density,  described  as: 


d()> 

a - h 

dx. 

i 


dDc 

_ J_ 

dx 

J 


(2.30) 


Here  a  is  the  electrical  conductivity  and  F  is  the  Faraday  constant.  The  electrical  conductivity  is 
calculated  as: 

N 

^  =  (2.31) 

y=i 

It  is  noted  that  in  the  solution  methodology,  the  electroneutrality  is  not  explicitly  applied. 
Instead,  transport  equation  is  solved  for  each  species  independently,  using  the  electrical  drift 
term  obtained  by  solving  current  conservation,  which  is  obtained  using  electroneutrality 
assumption. 


Approach  2  (Electroneutrality  not  Assumed) 

In  this  approach,  the  electrical  potential  is  obtained  by  solving  Poisson  equation: 

d  d(j) 
dx.  dx. 


z.c. 

j  j 


y=i 


(2.32) 


where  s  is  the  dielectrical  constant  of  the  medium.  With  this  formulation,  when  the  source  terms 
on  the  right  side  dominate,  problems  appear  in  obtaining  the  converged  solution.  This  issue  will 
be  addressed  in  future. 
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Model  validation 


In  order  to  verify  the  eurrent  numerical  approach,  the  code  is  validated  with  experimental  data  of 
Bharadwaj  (private  communication).  The  schematic  of  the  system  is  shown  in  Figure  2.13.  The 
system  used  consists  of  an  elliptical  shaped  geometry.  The  sample  is  prepared  in  low 
conductivity  buffer  (shown  in  gray)  and  the  high  conductivity  buffer  is  then  electrokinetically 
driven  from  the  buffer  well,  on  left,  towards  the  sample  well,  on  the  right  and  the  sample  stacks 
at  the  interface  between  the  two  electrolytes.  It  is  anticipated  that  the  sample  will  stack  at  the 
interface  between  the  two  electrolytes.  Three  species  can  be  used  in  this  model  problem  and  the 
parameters  used  are  listed  below  in  Table  2.2. 


Table  2.2:  Parameters  used  in  single-step  stacking  simulation. 


Species  name 

Diffusivity 

Mobility 

Charge 

A+ 

lE-9 

5.38E-13 

1 

B- 

7E-10 

3.174E-12 

-1 

C- 

5E-10 

3.12E-13 

-2 

8  cm 


^50  pni 


High 

buffer 


E  =  50  V/cm 

- ► 


T 


20  pm 


conductivity 


Low  conductivity 
sample 


Figure  2.13;  Schematic  of  single-step  stacking  device. 


The  hydrodynamic  properties  of  the  electrolyte  are  those  of  water.  The  EOF  mobility  is  set  to 
7E-10  m^/(V  s).  The  initial  concentration  is  set  to: 

C^=Q.5C^^{y  +  \-{r-\)erf{a{x-s)) 

Cc=Cc^{^  +  erf{a{x-s)) 

where  y=4  is  the  ratio  of  buffer  concentration  in  the  background  high  conductivity  and  to  that  in 
the  low  conductivity  sample.  In  the  simulation  the  following  parameters  are  chosen: 

s  =  \E-4,a  =  5EA,C^^=\E-XCc^=\E-6 
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The  concentration  of  species  B-  is  determined  from  the  electroneutrality  condition.  A  50V/cm 
electrical  field  is  applied  along  the  channel  and  simulates  the  evolution  of  sample  transport. 


Experiment  Simulation 

Figure  2.14:  Comparison  of  experimental  observation  and  simulation  of  sample  stacking  with 

EOF  induced  dispersion. 

The  snapshots  are  taken  at  t=0,  278,  558  and  834ms.  The  channel  width  is  50pm. 


A  2D  analysis  that  compares  the  images  (Figure  2.14)  shows  a  good  qualitative  agreement.  The 
team  then  simulates  with  an  assumption  of  an  infinite  supply  of  the  sample  and  the  interface 
moves  across  the  horizontal  channel.  The  first  set  of  simulation  was  performed  using  ID 
analysis  code.  The  results  from  this  simulation,  shown  in  Figure  2.15,  clearly  demonstrate  the 
stacking  of  specie  C  at  the  interface.  The  peak  value  is  almost  linear  as  observed  before  by  the 
Stanford  Group.  The  following  is  observed: 
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The  stacking  ratio  of  simulated  results  is  considerably  higher  than  that  of  experiments; 


•  The  translation  of  the  sample  peak  is  slower  than  that  of  the  experiment. 


Column  Position  (m) 

Figure  2.15;  Comparison  of  simulated  (dash  curves)  and  experimental  results  (solid  curves)  of  the 
normalized  concentration  at  the  centerline  of  the  channel. 


This  discrepancy  may  possibly  arise  from: 

1 .  The  electrostatic  interaction  between  analyte  molecules,  neglected  in  the  present  simulation, 
may  reduce  peak  height. 

2.  The  translation  speed  of  the  peak  is  a  function  of  the  electroosmotic  mobility,  which  may 
exhibit  more  complex  dependency  on  local  ionic  strength  and  pH. 

From  extensive  simulations  it  is  observed  that  the: 

•  Net  neutrality  assumption  under-predicts  electrical  field  near  the  electrodes  and 
subsequently  reduced  sample  stacking. 

•  Non-net  neutrality  assumption  predicts  a  higher  electrical  field  strength  leading  to 
enhanced  stacking. 

To  address  above  issues,  the  team  developed  models  for  concentrated  ionic  solutions  and 
variable  zeta  potential  models. 

(a)  Model  for  Concentrated  Ionic  Solutions: 

During  sample  stacking,  it  is  normal  to  stack  up  analytes  as  high  as  1000  times  their  initial 
concentration.  At  that  instant,  the  solution  becomes  highly  concentrated  and  analyte  particle- 
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particle  interaction  becomes  significant.  This  results  in  deviations  in  the  values  of  electrical  and 
physical  properties  from  their  respective  values  when  the  solution  is  dilute.  To  address  this 
issue,  the  team  has  developed  “unit-cell”  based  highly  concentrated  solution  models  to  simulate 
sample  stacking. 

(b)  Model  for  Variation  ofZeta  Potential  With  Respect  to  Ionic  Concentration: 

The  team  implemented  a  model  that  describes  variation  of  zeta  potential  as  a  function  of  ionic 
concentration  based  on  the  work  by  Scales,  et.  al.  ("Electrokinetics  of  the  Silica-Solution 
Interface:  A  Flat  Plate  Streaming  Potential  Study.  "The  ACS  Journal  of  Langmuir  Surfaces  and 
Colloids  8(3):  965-974).  The  accumulation  of  oppositely  charged  species  on  surface  will 
neutralize  the  surface  charge  it  carries  and  hence  reduce  the  zeta  potential  Such  a  coupled 
model  will  exhibit  different  behavior  from  non-coupled  model,  where  zeta  potential  is  assumed 
to  be  constant.  The  non-uniformity  of  zeta  potential  will  introduce  the  internal  pressure  gradient 
in  a  fashion  similar  to  that  seen  in  field  amplified  sample  stacking.  Such  internal  pressure 
gradient  will  induce  Taylor  dispersion  that  may  reduce  stacking  efficiency. 

(c)  Numerical  Simulation  of  Full  Poisson  Equation: 

Numerical  simulation  of  fully  coupled  Poisson  equation  is  extremely  challenging  when  there  is  a 
large  charge  accumulation  leading  to  very  large  source  terms  and  causing  a  numerical  stiffness 
problem.  In  this  regard  the  team  has  implemented  a  full  Poisson  equation  solver  that  uses  a 
SUPER  LU  scheme  to  solve  stiff  equation  (due  to  high  source  term  that  rises  from  charge 
accumulation)  that  describes  electric  field  distribution. 

2.1. 7  Model  for  EK  Phenomena  in  Concentrated  Electrolyte  Solutions 

During  several  electrochemical  processes  prevalent  in  the  biotech  industry  (e.g.,  sample 
stacking)  sample  concentration  increases  and  the  solution  may  become  highly  concentrated. 
When  this  happens,  the  dilute  approximation,  based  on  which  current  ACE-i-  models  are 
formulated,  may  become  unstable/invalid.  Hence,  the  team  has  developed  methodologies  for 
solving  highly  concentrated  solution  based  on  a  semi-implicit  coupled  approach.  This  is 
discussed  next. 

Stability  of  the  concentration  and  potential  solutions  is  first  sought  to  be  improved  by 
strengthening  the  coupling  during  numerical  solution.  CFD-ACE-i-  solves  charged  analyte 
transport  equations  under  electrostatic  field  as  following: 

VfiV  ^  =  F'^  zf. 

5c  V  V 

^  +  V  •  (VC,)  +  V  •  {z,co,Ecf  =  V  •  (D,Vc,)  (2.33) 

ot 

E  =  -V(l> 
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The  original  implementation  will  face  convergence  problem  if  the  electro  non-neutrality  term 
dominates.  The  previous  equations  are  manipulated  multiplying  the  valance  and  Faraday 
constant  and  summing  through  all  analytes: 


^^fSfiV*  +  /r^2V.(yc,)  +  V.(Fyz;-«,ic,)  =  FX--V'(AVc,>  (2.34) 

This  equation  gives  the  charge  density  at  next  time  step  in  terms  of  ion  concentration  at  old  time 
step: 

(FX^)C,r=(FX^)C,r-V•(A^FXz>,.^c,)  +  A^F{X^)[V•(AVc,)-V•(vc,)]}(2.35) 
The  modified  electrostatic  equation  for  new  time  step  is  therefore, 

V  ■  («  +  A<  ■  Fy  )  V«t  =  (Fy  z,c, )”  +  A(  ■  {y  z,  [  V .  (Z),  Vc,)  -  V .  (VC,)])  (2.36) 

For  steady  state  simulations,  the  time  step  can  be  simply  set  a  large  number.  This  methodology 
has  been  developed  (by  CFDRC  under  a  different  contract)  and  has  been  used  in  this  project  to 
simulate  FASS.  The  numerical  scheme  is  summarized  in  Figure  2. 16. 
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Figure  2.16:  A  semi-implicit  scheme  for  Poisson  equation  coupled  with  species  transport. 


2.1.8  Microsphere-Based  Sandwich  Immunoassay  Modeling  Capability 

The  basic  concept  of  tagged  antibodies  (primary  antibody)  anchored  to  the  microbead  surface 
and  then  exposed  to  target  analytes  in  the  flow  system  is  well  known.  In  order  to  quantify  the 
antigen  bound  and  also  to  increase  the  sensitivity  of  the  signals  and  accuracy  of  measurement, 
“sandwich  immunoassays”  is  used.  Here,  one  more  step  is  added  to  the  immunoassay,  which  is 
known  as  the  “Reporter  Step”.  These  types  of  assays  have  found  wide  application  in 
biotechnology.  Figure  2.17  below  explains  the  basic  concept  involved  in  sandwich 
immunoassays. 
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Figure  2.17:  Schematic  of  the  sandwich  immunoassay. 

(a)  Primary  antibody  coated  microsphere,  (b)  antigen  bound  to  primary  antibody  on  microsphere, 
and  (c)  reporter  step. 


In  this  case  the  primary  antibody  is  immobilized  on  the  microsphere.  Buffer  containing  antigen 
is  added  and  allowed  to  react  with  the  bound  primary  antibody  as  shown  in  Figure  2.17b.  After 
sufficient  mixing  and  binding,  excess  buffer  and  antigen  are  washed  off  After  this  step,  the 
secondary  antibody  (specific  for  the  distinct  epitope  on  the  antigen)  is  added  and  allowed  to  react 
with  the  bound  antigen.  Again,  after  sufficient  mixing  and  binding,  the  excess  secondary 
antibody  sample  is  washed  off  and  the  beads  are  sent  to  a  cytometer,  where  the  bead 
fluorescence  (on  excitation  of  laser)  is  measured.  The  intensity  of  fluorescence  on  the 
microsphere  surface  is  proportional  to  the  concentration  of  secondary  antibody  in  each  sample. 

Model  Validation  for  Sandwich  Assay 

The  sandwich  immunoassay  is  easily  modeled  using  the  existing  generalized  surface  chemistry 
solver  implemented  in  the  CFD-ACE+  multiphysics  engine.  During  a  simulation,  difftisional 
mass  transfer  to  the  bead  surface,  and  the  surface  reactions  are  balanced  for  all  surface  reactions 
simultaneously.  The  model  was  validated  first  against  semi-analytic  expressions  possible  for  a 
static  incubation  assay.  Here  Primary  Antibody  (PAb)  coated  microspheres  are  incubated  (no 
flow)  with  excess  of  Antigens  1  &  2  (100  nM).  The  set  of  reactions  are: 

Agl  +  FAb^Agl-FAb 

Ag2  +  FAb  ^  Agl  -  FAb 

SAb  +  Agl  -  FAb  «  SAb  -  Agl  -  FAb 

Model  equations  reduce  to  a  system  of  linear  ODEs  for  this  static  incubation  case,  which  were 
solved  with  using  the  ordinary  differential  equation  solver  (ode45)  available  in  MATLAB®  and 
the  results  are  compared  with  CFD-ACE-i-  results  for  the  static  case.  The  time  evolutions  of  all 
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three  surface  species  are  identical  (Figure  2.18),  validating  the  integration  of  multi-step  kinetic 
models. 


(a)  (b) 

Figure  2.18:  (a)  Static  incubation  assay  with  microsphere. 

(b)  ACE+  binding  results  compared  with  the  numerical  solution  of  the  underlying  ODE’s. 


Demonstration  of  the  Microsphere-Based  Immunoassay  Model 

A  demonstration  of  the  multi-step  bead  assay  model  in  a  typical  lab-on-a-chip  as  illustrated  in 
Figure  2.19.  Primary  antibody  coated  beads  are  injected  from  Inlet  1.  A  mixture  of  antigens  is 
injected  from  Inlet  2.  The  antigen  mixes  with  antibody  coated  antibodies  in  the  vertical  charmel. 
The  excess  antigens  are  washed  out  via  Outlet  1 .  Next,  fluorescence  tagged  secondary  antibodies 
are  injected  via  Inlet  3.  The  secondary  antibody  complexes  only  with  the  antigen  that  is  already 
captured  on  the  beads.  Excess  secondary  antibody  is  washed  out.  Finally,  the  amount  of 
fluorescence  detected  on  the  beads  (via  the  bound  secondary  antibody)  serves  to  quantify  the 
presence  of  the  original  antigen. 
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Antigen  coaled  hegils 


Figure  2.19;  Schematic  of  multi-step  microsphere  assay  demonstration/validation  study 


A  demonstration  assay  was  carried  out  using  values  based  on  the  Interleukin  -  receptor  system. 
Figure  2.20a  shows  a  snapshot  of  beads,  covered  with  the  primary  antibody  entering  the  turn 
(leading  to  Outlet  1).  While  a  small  amount  of  beads  are  lost  to  the  outlet,  the  vast  majority  make 
it  into  the  continuing  channel.  Peak  coverage  of  about  40%  is  seen  indicating  that  the  length, 
while  sufficient,  is  not  enough  to  saturate  the  beads.  The  secondary  binding  is  shown  in  Figure 
2.20b.  Once  again,  very  few  beads  are  lost  in  the  turn.  Interestingly,  less  than  a  third  of  the 
primary  covered  antibody  is  bound  to  the  secondary,  signal-generating,  antibody  indicating  that 
this  particular  assay  can  indeed  be  further  optimized  using  the  existing  simulation  models. 


Figure  2.20:  (a)  Primary  antibody  binding,  (b)  Secondary  antibody  binding. 
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Demonstration  of  a  Multiplexed  Bead  Assay  Application 


The  multiplexed  bead  assay  is  demonstrated  on  a  T-junction  (commonly  found  in  many 
microfluidic  lab-on-a-chip  systems)  to  mix  microbead  solutions  with  sample  solutions  (Figure 
2.21).  The  horizontal  channel  was  100  pm  wide  and  2200pm  long,  while  the  vertical  channel 
was  100  pm  wide  and  500  pm  long.  Two  different  analytes  of  the  cytokine  family,  IL-2  and  IL- 
4  were  used  for  these  simulations.  The  diffusivity  of  both  the  analytes  was  set  at  lO  ’*  m^/s.  10 
pm  beads  coated  with  the  respective  antibody  against  IL-2  and  IL-4  were  used.  The  adsorption 
constant  (Ka)  for  the  interaction  between  IL-2  and  its  antibody  was  8x10^  M'*  s'*  and  between 
IL-4  with  its  corresponding  antibody  was  5x10  M'  s'  .  The  total  available  site  for  binding  was 
10'^  moles/m^.  IL-2  and  IL-4  plugs  were  injected  for  3  seconds  each,  in  the  horizontal  channel  at 
a  velocity  of  0.001  m/sec  with  the  concentration  of  lO'*’  M.  IL-2  was  injected  first  followed  by 
IL-4  with  an  interval  of  6  seconds  between  the  injections.  The  beads  coated  with  the  respective 
antibodies  were  injected  from  the  vertical  channel  in  a  buffer  solution  at  the  same  velocity.  The 
beads  were  injected  3.5  seconds  after  the  introduction  of  IL-2,  timed  so  that  both  beads  and  the 
cytokines  arrive  at  the  junction  at  the  same  time.  The  total  time  of  the  assay  was  30  seconds. 


Figure  2.21  shows  the  channel  at  5  seconds  with  the  first  analyte  (IL-2)  and  bead  plugs  (coated 
with  antibodies  to  both  IL-2  and  IL-4)  entering  the  respective  channels.  As  expected  the  IL-2  is 
diffusing  across  the  channel  and  the  beads  show  no  coverage  as  they  have  not  yet  encountered 
the  cytokine.  At  15  seconds  (Figure  2.22),  IL-4  has  entered  the  channel,  while  the  beads  are 
mixing  with  the  IL-2  plug.  Note  that  only  some  of  the  beads  in  contact  with  the  IL-2  are 
changing  color,  indicating  the  presence  of  IL-2  specific  antibody  while  the  IL-4  antibody  coated 
beads  remain  colored  blue. 


Figure  2.23  shows  the  channel  at  18  seconds  when  IL-2  is  exiting  the  channel  and  IL-4  has 
reached  the  center  of  the  channel.  Notice  that  some  of  the  beads  (now  the  ones  coated  with  IL-4 
specific  antibody)  are  detecting  the  presence  of  IL-4  while  no  change  is  seen  in  the  IL-2  antibody 
coated  beads.  At  a  much  later  time  (Figure  2.24)  IL-2  has  exited  the  channel  and  the  IL-4 
antibody  coated  beads  are  seen  to  have  detected  the  presence  of  IL-4.  This  study  can  be  readily 
extrapolated  to  include  an  arbitrary  number  of  analytes  and  beads. 
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Figure  2.21:  T  junction  at  time  =15  seconds  following  injection. 


Figure  2.22:  T  junction  at  time  =15  seconds  following  injection. 


Figure  2.23:  T  junction  at  time  =18  seconds  following  injection 
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Beads  interacting  with  IL-4 
(Blue  indicate  beads  with  antibody  to  IL-2  and 
the  rest  indicated  beads  with  antibody  to  IL-4, 
red  being  maximum  coverage) 


Figure  2.24:  T  junction  at  time  =  25  seconds  following  injection  (IL-2  has  been  washed  out  of  the  channel 

at  this  time) 

Effect  of  Analyte/Bead  Flow  Rates  on  Detection 


The  goal  of  all  biochemical  assays  is  to  yield  the  maximum  signal  in  the  shortest  time  possible. 
The  signal,  however,  depend  upon  a  multitude  of  variables  such  as  geometry,  flow  rates, 
concentrations  etc.  in  a  very  complex  and  intertwined  manner.  Here,  the  team  has  chosen  to 
study  the  effect  of  one  of  the  most  important  and  easily  controllable  parameters  i.e.  antigen  and 
bead  flow  rates.  The  geometry  considered  is  a  Y-junction  employed  in  many  microfluidic 
systems.  The  length  of  the  arm  of  channel  is  0.5  mm  and  its  width  is  0.05  mm.  The  length  of  the 
straight  channel  is  2  mm  and  its  width  is  0.1  mm.  The  sample  analyte  (plug)  and  buffer  enters 
from  upper  arm  of  the  channel  and  beads  in  a  buffer  solution  enter  from  lower  arm  of  channel. 
The  diffusivity  of  the  analyte  is  set  at  5x10'  m  /s.  The  bead  fluorescence  (coverage)  is 
measured  at  the  exit  of  the  channel  (at  outlet).  The  analyte  concentration  used  in  all  the  cases  is 
10'*^  M.  5  pm  beads,  coated  with  receptors,  are  released  in  the  lower  arm.  An  adsorption  constant 
of  10^  is  prescribed,  for  the  purposes  of  this  calculation. 


Figure  2.25:  Analyte  concentrations  and  bead  coverages  for  a  flow  ratio  of  5(An):5(B). 
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For  the  baseline  case,  equal  flow  rates  of  analyte  and  bead  entering  the  channel  are  considered 
(Figure  2.25).  Note  that  the  analyte  and  bead  each  occupy  half  the  channel  with  very  little 
mixing.  The  suboptimal  mixing  and  presentation  of  sample  to  the  bead,  causes  a  low  level  of 
surface  coverage.  Considering  that  the  threshold  of  detection  in  many  cases  requires  around  30% 
coverage  and  the  present  protocol  shows  substantially  less  coverage  than  that,  it  can  safely  be 
concluded  that  the  test  as  presently  conducted  has  a  high  probability  of  yielding  a  false  negative. 


Figure  2.26;  Analyte  concentrations  and  bead  coverage  for  a  flow  ratio  of  5(An);l(B) 


Decreasing  the  bead  flow  rate  can  be  intuitively  expected  to  increase  the  residence/contact  times 
between  the  sample  and  antibody  (on  the  surface  of  the  beads)  and  thereby  increase 
coverage/detection.  This  is  indeed  seen  to  happen  as  seen  in  Figure  2.26.  Moreover,  the 
hydrodynamic  effect  of  focusing  the  beads  in  the  low-speed  fluid  near  the  wall,  also  enhances  the 
signal  via  increased  diffusive  transport  of  analyte  to  the  surface  of  the  beads 


The  overall  signal  evolution  with  time  at  the  detection  point  (at  the  exit  of  the  channel)  is  shown 
in  Figure  2.27.  Firstly,  as  expected,  the  time  to  detect  is  shortest  in  the  5:5  case,  as  the  beads 
travel  fastest  through  the  system  and  reach  the  detection  point.  However,  as  seen  earlier,  the 
coverage  levels  are  low  with  the  possibilities  of  false  negatives,  high.  The  5:1  case,  where  the 
beads  travel  slowly,  on  the  other  hand  exhibits  a  higher  coverage  that  is  also  sustained  for  longer 
times.  Hence  it  is  concluded  that  flowing  the  beads  in  slowly  is  the  preferred  mode  of  operation 
in  this  assay. 
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Figure  2.27:  Signal  (sensorgrams)  evolution  with  time. 


2.1.9  Simulation  of  Combined  DEP  and  Microbead  Assay 


The  DEP  and  the  microbead  assay  models  were  successfully  coupled  to  simulate  analyte  binding 
on  a  microbead  cluster  while  they  were  being  transported  across  a  microchannel  to  due  AC 
electric  field.  For  this  simulation,  fluid  flow  is  solved,  AC  electric  field,  DEP-induced  bead 
motion  and  analyte-antibody  binding. 

A  total  of  100  microspheres  were  introduced  randomly  at  the  intersection  of  a  cross-channel 
configuration  that  is  100  x  100  pm  in  dimension.  Water  is  considered  as  the  fluid  medium 
occupying  the  cross  channel  initially.  The  physical  and  electrical  properties  of  the  medium  and 
spheres,  provided  by  UC  Santa  Barbara,  are  listed  in  table  2.3. 

Table  2.3:  Physical  and  electrical  properties  of  the  medium  (DI  water)  and  microspheres  (by 


UCS] 

B) 

Frequency  (Hz) 

100  Hz 

Fluid  Density 

1000  Kg W 

Particle  Radius 

0.1  pm 

Particle  Density 

1000  KgW 

Conductivity  of  Particle 

l.OE-14  1/Ohm-m 

Conductivity  of  the  Fluid 

0.56  1/Ohm-m 

Relative  Permittivity  of 
Particle 

2.61 

Relative  Permittivity  of 
Fluid 

80 
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Microspheres  were  coated  with  IL-2  antibody.  IL-2  analyte  was  assumed  to  be  initially 
restricted  to  the  left  and  right  chamber  of  the  domain.  The  diffusivity  of  IL-2  was  set  at  l.OE-10 
m  /s.  The  adsorption  constant  (Ka)  for  the  interaction  between  IL-2  and  its  antibody  was  8x10 
M  '  s  '.  The  available  binding  site  density  on  the  microspheres  was  lE-06  moles/m^.  An  AC 
electric  field  of  10  V  is  applied  to  the  top  and  bottom  electrodes.  The  frequency  of  the  AC  field 
was  kept  at  100  Hz.  Since  the  relative  permittivity  of  the  particle  is  less  than  that  of  the  medium, 
it  will  experience  a  negative  DEP  (moves  toward  electric  field  minima). 

Figures  2.28  and  2.29  show  the  initial  concentration  contours  of  IL2,  and  the  contours  of  electric 
field  squared,  respectively.  The  random  locations  of  the  microsphere  are  superimposed  the  top 
of  the  contour  plot.  Figure  2.29  illustrates  the  sequence  of  operation  as  it  occurs  in  the  cross 
channel.  From  Figures  2.29a  to  2.29d,  the  beads  migrate  horizontally  on  either  direction 
(positive  and  negative  x-direction)  towards  the  electric  field  minima  (Figure  2.29a)  as  expected. 
As  they  migrate,  they  displace  fluid  volume  in  the  left/right  chambers,  which  sets  up  chaotic 
mixing  regions  within  the  domain.  The  IL-2  in  the  horizontal  channel  also  starts  to  bind  with  the 
IL-2  antibody  coated  on  the  beads.  This  is  illustrated  via  a  color  change  on  the  bead  (blue  to 
red).  At  time  t  =  2  seconds,  almost  all  the  beads  have  accumulated  near  the  electric  field  minima 
and  the  antigen-antibody  binding  is  almost  complete. 

IL2  L2E+010 

9.9E-007_ 

9.0E-007 
3-0E-0O7 
TOE -007 
6.0E-007 
5.0E-0O7 
4.DE-D07 
:5  0E-007 
2.0E-007 
l.DE-007 


S.4E-Oia  9-1E  +  004 


Figure  2.28a:  Concentration  contours  of  IL2  at  Figure  2.28b:  Contours  of  electric  field  squared  at 
time  t  =  0  time  t  =  0 
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(a)  t  =  0.5  sec 


(b)  t  =  1.0  sec 


(c)t=1.5sec  (d)t  =  2.0sec 

Figure  2.29;  Location  of  microspheres  shown  at  different  time  levels. 


The  color  of  the  sphere  indicates  extend  of  antigen-antibody  binding  (blue  ->  all  receptor  sites 
available,  red  ->  no  site  available).  Contour  levels  of  1L2  are  also  shown. 


2.1.10  Architecture  for  Biomolecular  Databases 


Primary  integration  of  the  Biochemistry  and  Electrochemistry  modules  in  CFD-ACE+  is 
complete.  The  Electrochemistry  module  has  a  detailed  chemistry  database  structure  set-up  for 
gas-phase  chemistry  for  semiconductor  applications.  The  team  designed  similar  database 
architecture  for  the  biomolecular  reaction  kinetics  data  as  well.  A  sample  graphical-user- 
interface  for  specifying/choosing  biomolecular  surface  interactions  is  shown  below  in  Figure 
2.30.  After  completion  of  this  task,  the  user  can  readily  choose  from  a  wide  range  of  data 
pertaining  to  biomolecular  kinetics,  presented  in  an  intuitive,  easy-to-use  format. 
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Figure  2.30:  Developed  architecture  for  surface  biokinetics  manager. 

This  user  interface  keeps  track  of  reactions,  reagents,  and  the  salient  reaction  kinetics  and 

transport  parameters. 

2.1.11  Neural  Net  Modeling  for  Data  Extraction 


Work  was  performed  on  the  adaptation  of  the  Artificial  Neural  Network  (ANN)  models  to  the 
extraction  of  the  coefficients  associated  with  surface  reaction  kinetics.  The  approach  taken  is 
based  on  treating  the  ANN  as  a  reduced  model  for  time  series  prediction.  The  time  series 
represents  the  time  history  of  surface  concentration  as  a  function  of  input  flaw  rate,  temperature, 
as  well  as  the  maximum  allowable  surface  and  initial  bulk  concentrations.  System  geometry  is 
taken  as  constant.  This  non-linear  dynamical  system  is  being  represented  by  a  static  topology 
ANN  in  a  dynamic  configuration  of  the  general  form: 


y(t+l)  =  f[y(t),  y(t-l), ...,  y(t-m),  u(t),  u(t-l), ...,  u(t-n)] 

(2.37) 


where  y(t),  u(t)  represent  the  time  dependent  output  and  input,  respectively.  The  work  was 
performed  in  three  stages:  (1)  General  ANN  preparation;  (2)  Input/Output  Data  Preparation;  and 
(3)  Model-Based  ANN  preparation. 
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1)  General  ANN  Preparation:  The  static  ANN  is  cast  in  a  general  feed- forward  Multi  Layer 
Perceptron  (MLP)  topology  with  time  lagged  elements  and  global  feedback  capability.  The  work 
involved  using  model  problems  for  testing  of  training  techniques  for  selected  candidate  ANN 
topologies.  The  topologies  under  consideration  include  ANNs  with  feedback  elements.  Such 
representations  can  in  principle  model  very  complex  dynamical  systems  but  require  more  general 
training  techniques.  Test  of  these  ANNS  are  nearly  complete.  The  final  selection  of  the  most 
suitable  ANN  architecture  will  be  based  on  test  results  using  experimental  data. 


2)  Input/Output  Data  Preparation:  Since  the  time-dependent  process  that  the  ANN  is  modeling  is 
of  finite  duration  all  ANN  training  will  be  done  off-line  using  archived  experimental  and 
computational  data.  The  existing  ANN  software,  designed  for  on-line  computations  using  real¬ 
time  input,  was  extended  to  accept  archived,  file-based  data.  Tests  of  this  extension  are  complete. 


3)  Model-Based  ANN  Preparation:  A  general  model-based  ANN  implementation  for  extraction 
of  the  reaction  rate  coefficients  has  been  completed  and  tested  on  simple  problems.  The 
implementation  involved:  (a)  modification  of  the  ANN  structure  to  allow  specification  of 
analytical  expressions  to  compute  the  final  ANN  output  (concentrations);  and  (b)  modification  of 
the  training  procedure  to  include  the  effects  of  model  coefficients  variation  on  the  output  error. 
Complete  tests  and  training  of  this  formulation  will  be  performed  using  the  topology  as  described 
in  the  paragraph  above. 


ANN  Training  Set  Generation 

Training  data  for  the  neural  network  was  generated  using  a  model  of  the  SPR  biosensor  - 
Biacore.  The  data  sets  were  generated  using  the  CFD-ACE-i-  multiphysics  software  to  simulate 
the  effects  of  fluid  flow,  mass  transport  via  diffusion,  and  the  binding  of  analyte  to  receptor  on 
the  lower  surface  of  the  Biacore  device.  The  surface  binding  module  assumes  isothermal, 
Langmuir  adsorption  through  an  adsorption  (Ka)  and  desorption  (Kd)  term.  The  system  variables 
held  constant  were  the  sensor  surface  receptor  concentration  of  2.5  E-8  mole/m^,  inlet  analyte 
concentration  of  225  nM,  flow  rate  of  50  microL/min.  The  analyte  was  flowed  in  for  50  seconds 
followed  by  a  30  second  wash  phase  to  mimic  actual  experimental  conditions.  Kinetic  constants 
Ka  and  Kd  were  paired  over  a  wide  range  to  capture  the  effects  of  kinetic  and  mass  transport 
limitations  (Figure  2.31),  so  as  to  assure  the  validity  of  ANN  results  over  the  widest  possible 
range  of  experimental  conditions.  100  data  sets  were  generated  for  the  neural  network  training. 
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Data  Set  Parameter  Space 


Figure  2.3 1 :  Parameter  space  for  training  ANN 
(Kinetic  to  diffusion  controlled.  Adsorption  to  desorption  dominant) 

Neural  Network  Training 

Artificial  Neural  Networks  (ANN)  with  the  feedforward  Multi  Layer  Perceptron  (MLP) 
architecture  are  being  used  in  this  study.  The  ANNs  are  trained  to  correlate  the  time  variation  of 
surface  concentration  C(t)  to  the  numerical  values  of  the  kinetic  coefficients  (Ka,  Kd),  over  a 
range  of  the  latter.  The  networks  thus  define  a  reduced  model  for  a  physical  process,  capable  of 
predicting  the  time  variation  of  C(t)  for  selected  (Ka,Kd)  values.  Currently  the  nets  are  being 
trained  using  the  standard  Backpropagation  (BP)  algorithm  with  a  momentum  correction. 
Typical  values  of  the  learning  rate  and  the  momentum  factor  are  on  the  order  1.  Several 
multiplayer  topologies  are  being  investigated  with  inputs  and  outputs  defined  by  (t,Ka,Kd)  and 
C(t),  respectively.  The  training  data  contains  8000  datapoints  representing  100  concentration 
histories  for  different  (Ka,  Kd)  values,  with  80  time  steps  per  history.  Tests  to  date  indicate  that 
the  following  topologies  give  favorable  results:  3-20-10-1  and  3-20-10-10-1  The  first  and  last 
numbers  in  the  previous  strings  represent  the  number  of  inputs  and  outputs,  respectively,  and  the 
intermediate  values  denote  the  number  of  neurons  in  the  hidden  layers. 
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Parameter  Extraction  Results  for  Neural  Network 

A  fully  trained  neural  network  allows  for  rapid  extraction  of  kinetic  constants  using  a  non-linear 
least  squares  optimization  protocol.  The  experimental  data  is  essentially  curve  fit  using  the 
neural  network  as  the  generator  of  the  theoretical  curve.  The  neural  network  takes  as  input  a 
(Ka,  Kd)  pair  and  a  set  of  discrete  time  points  and  gives  as  output  the  surface  concentration  at 
each  time  point.  In  the  current  case,  a  Gauss-Newton  non-linear  least  squares  algorithm  is  used 
to  minimize  the  squared  error  between  the  experimental  data  and  the  neural  network  generated 
data.  The  “experimental  data”  used  in  this  case  was  a  data  set  was  generated  with  CFD-ACE+ 
for  a  Ka  =  8E+6  l/(M-s)  and  Kd  =  0.24  1/s.  This  pair  was  not  one  of  the  pairs  in  the  original 
training  set,  but  was  contained  within  the  bounds  of  the  training  set. 

The  kinetic  constants  extracted  from  the  experimental  data  were  Ka  =  7.7E+6  l/(M-s)  and  Kd  = 
0.19  1/s,  agreeing  well  with  the  known  parameters  (Figure  2.32).  The  time  for  extracting  the 
kinetic  parameters  from  the  data  sets  was  approximately  5  minutes  using  the  neural  network. 
Please  note  that  parameter  extraction  using  a  fully  three  dimensional  ACE+  simulations  takes 
approximately  15-30  hours  depending  on  the  initial  guess.  The  time  for  training  data  set 
generation  was  approximately  6  hours,  with  a  training  time  of  1-3  hours  depending  on  the 
training  algorithm. 
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ACE+  «  2  hours  (2D) 
NN  5  minutes 


Figure  2.32:  Comparison  of  ACE+  and  ANN  rate  constants,  along  with  the  time  taken. 


(Training  time  is  around  1-3  hours). 


2.1.12  Proteomics  -  Array  Based  Biosensors  Models  and  Validation  (CFDRC  and  Univ.  of 
Utah) 


In  many  of  the  SPR  sensors  (including  Biacore)  the  immobilized  species  (receptor)  is  covalently 
attached  to  a  very  thin  hydrogel  matrix  (usually  a  dextran  matrix)  relative  to  the  height  of  the 
microchannel.  For  example,  in  the  Biacore  flow  cell  the  height  of  the  flow  channel  is  usually  on 
the  order  of  50  pm,  while  the  height  of  the  hydrogel  layer  is  generally  on  the  order  of  100  nm. 
This  scale  disparity  precludes  any  explicit  simulation  of  the  gel  layer.  However,  the  presence  of 
these  capture  surfaces  has  been  shown  to  be  of  critical  importance  and  must  be  taken  into 
account  in  high-fidelity  computational  models[30]. 
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The  most  common  method  used  for  modeling  ligand  transport  within  a  porous  layer  is  the 
cylindrical  pore  model  originally  proposed  by  [31].  This  model  accounts  for  the  equilibrium 
partitioning  of  the  bulk  concentration  into  the  hydrogel  layer  along  with  the  hindered  diffusion 
experienced  by  a  ligand  molecule  traveling  through  a  network  of  cylindrical  pores.  These 
partition  and  hindered  diffusion  coefficients,  in  turn,  depend  on  hydrogel  and  the  protein 
molecular  characteristics.  Based  on  prior  investigations,  [30],  [32]  the  team  have  developed  and 
implemented  within  the  convective-diffusive-reactive  framework  of  CFD-ACE-i-,  a  model  to 
explicitly  account  for  the  partitioning,  hindered  diffusion  and  volumetric  binding  effects  due  to 
the  gel  layer. 

The  main  difference  between  this  approach  and  prior  surface  biochemical  models  is  that  the 
receptor,  analyte,  and  receptor/analyte  complex  are  described  volumetrically  within  the  dextran 
matrix  (as  opposed  to  surface  concentrations).  Complex  effects  such  as  change  in  pore  size  due 
to  analyte-receptor  complexation  are  accounted  for  in  this  approach.  Partitioning  of  the  analyte 
from  the  bulk  solution  into  the  dextran  matrix  is  accomplished  following  the  vertical  cylindrical 
pore  model  of  Giddings  et  al.  (1968).  Here,  as  a  starting  point,  it  is  assumed  that  the  partitioning 
can  be  characterized  using  geometric  considerations  without  introducing  thermodynamic 
parameters.  The  dynamic  porosity  of  the  dextran  matrix  can  be  computed  as  follows: 


s{t)  =  1 


dextran  J 
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(2.38) 


Where  Na  is  Avagadro’s  Number,  bdextran  is  the  thickness  of  the  dextran  layer,  CVeceptor  is  the 
available  number  of  receptor  sites,  Ccompiex(t)  is  the  time  varying  analyte/receptor  complex,  and 
fAnaiyte  and  rReceptor  are  the  analyte  and  receptor  hydrodynamic  radii  respectively.  Now,  the 
average  vertical  pore  radius  can  be  estimated  as: 


'  Pore 


yl  -  S{t) 


3  3 

T  —  T 

Re  ceptor  Analyte 


C. 


Complex 


^  Re  ceptor  J 


'  Re  ceptor  '  Analyte 


C. 


Complex 


w 


r° 

^Receptor  J  J 


The  partitioning  coefficient  is  then  specified  as: 
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The  concentration  of  unbound  analyte  within  the  dextran  matrix  (no  gradients  are  assumed  to 
exist  internal  to  the  dextran  layer). 
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Cdextran  fj\  Interface  ( j\i  ( 

Analyte\t)  =  C  Analyte  VWPartitionVm) 


(2.41) 


The  boundary  conditions  in  the  CFD  model  need  to  be  modified  to  include  the  presence  of  the 
dextran  layer  and  the  partitioning  that  takes  place.  Figure  2.33  illustrates  the  adjusted  flux 
balance  boundary  condition  that  accounts  for  the  presence  of  the  dextran  layer.  The  final 
balancing  of  the  transport  and  reactive  fluxes  is  of  the  form: 

- - - (cf'^  ^  =0  (2.42) 

^  ^  Analyte  Analyte  }  a  Complex  a  \  IT  K I  Analyte  Acceptor 

^  ^dextran 

The  first  term  refers  to  the  diffusive  flux  of  analyte  between  the  cell  center  and  the  dextran 
interface.  The  second  and  third  terms  describe  desorption  and  adsorption  of  the  analyte/receptor 
complex,  respectively. 


Figure  2.33:  Illustration  of  the  adjusted  flux  balance  boundary  condition  accounting  for  the  presence  of 

the  dextran  layer. 

The  diffusive  flux  is  balanced  from  the  center  of  the  cell  to  the  dextran  interface.  The  interfacial 
concentration  is  then  assumed  to  partition  into  the  dextran  matrix. 

The  above  model  has  been  implemented  in  the  previously  available  biochemical  reaction 
framework  inside  CFD-ACE+.  The  inputs  to  the  model  (in  addition  to  the  traditional  inputs)  are 
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the  hydrodynamic  radii  for  the  analyte  and  receptor,  the  thickness  of  the  dextran  layer  and  initial 
pore  radii  and  porosities.  An  example  calculation  (to  gain  insight  into  model  behavior)  was 
carried  out  using  the  full  3D  Biacore  geometry  (created  previously).  The  buffer  flowrate  was  set 
at  a  typical  value  of  100  pL/min  with  an  analyte  inlet  concentration  of  0.20  pM.  The  binding 
kinetics  were  set  at  Ka=57000  (M-s)''  and  Kd=0.045  s’'.  The  density  of  available  binding  sites 
was  set  to  9.7E-8  mol/m^.  The  key  parameters,  in  addition  to  the  Damkohler  and  Kinetoc 
numbers  described  earlier  for  the  virgin  surface  capture  simulations  described  earlier,  are  the 
ratio  of  analyte  to  pore  radii  and  the  overall  porosity.  As  expected  the  capture/coverage  decreases 
as  the  ratio  of  analyte  size  to  pore  size  increases  (Figure  2.34).  This  is  mainly  due  to  the 
hindrance  effects  associated  with  the  decrease  in  the  average  vertical  cylindrical  pore  radius  as 
the  pores  fill  with  analyte-receptor  complex.  This  is  mechanistically  reflected  by  the  partitioning 
coefficient  decreasing  with  decrease  in  pore  size  (alternatively  increasing  analyte  size). 
Characterizations  of  other  parameters  such  as  porosity  are  currently  underway. 


Figure  2.34:  Change  in  the  coverage  curve  with  the  ratio  of  analyte  radius  (Ra)  to  initial  vertical 

cylindrical  pore  radius  (Rp). 

This  calculation  was  performed  using  a  constant  receptor  radius  of  5  nm  and  a  dextran  thickness 
of  100  nm.  The  calculation  was  performed  on  the  full  3D  Biacore  chip  geometry. 


Model  data  is  available  from  U.  Utah  with  various  surface  treatments.  In  their  experiments,  three 
different  surface  treatments  were  experimented  with  three  different  monoclonal  antibody 
systems.  Their  results  (Figure  2.35)  showed  that  the  behavior  of  the  monoclonal  antibody 
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Amine  Dextran  Capture  Dextran  Absorbed  Gold 


systems  is  quite  different  on  the  different  surfaces.  For  example,  Table  2.4  lists  the  fitted 
association  constants  for  the  individual  system.  Note  the  wide  variation  in  kinetic  constants  that 
has  been  observed 


mAb  M2 


mAb  M2 


mAb  M6 


mAb  M6 


mAb  M7 


0  100  200  300  400  500  600  700 


Figure  2.35:  Illustration  of  the  differing  rates  of  adsorption  on  various  monoclonal  antibody  systems  and 

surfaces. 

The  three  surfaces  are  bare  gold,  amine  dextran,  and  capture  dextran.  The  kinetic  coefficients  for 
each  of  the  above  graphs  are  given  below  in  Table  2.4.  Data  were  obtained  using  the  Biacore-SPR 
device. 
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Table  2.4:  Binding  rate  eonstants  determined  from  SPR  analysis  of  PSA/antibody  interactions 


Absorbed  Gold 
Cature  Dextran 
Amine  Dextran 


Association  Rate  (Ka) 
mAb  M2  mAb  M6  mAb  M7 

1.36E+06  1.36E+05  1.64E+06 

8.05E+05  1.07E+05  5.22E+05 

2.55E+05  6.51  E+04  2.42E+05 


2.1.13  Design  of  New  Biacore  Flow  Cell  (CFDRC  and  Univ.  of  Utah) 

CFDRC  worked  with  David  Myszka  in  the  design  of  a  new  flowcell  for  the  Biacore  instrument 
whose  purpose  is  to  examine  biomolecular  interactions  in  the  presence  of  mass  transport  limited 
conditions.  In  their  planned  design,  this  is  accomplished  using  a  recess  in  the  flow  cell  before 
and  after  the  sensor  patch.  The  team  is  carrying  out  full-fledged  sensing  simulations  to  assess  the 
impact  of  the  recess  on  flow,  analyte  transport  and  binding  within  the  new  flow  cell.  The  flow 
channel  domain  (with  the  notch)  is  shown  in  Figure  2.37.  Figure  2.37b  shows  the  notch  along 
with  associated  concentration  profiles.  The  flow  extends  inside  the  notch,  slowing  down  as  it  fills 
the  notch  and  also  driving  small  secondary  vortices.  Possible  implications  for  analyte  transport 
are  twofold  (a)  there  is  a  time-lag  or  slow  down  associated  with  the  longer  time  taken  to  traverse 
the  notch  and  (b)  dead  zones  in  the  notch  comers  may  pose  potential  lingering  cross 
contamination  problems. 


Figure  2.36:  (a)  Flow  cell  with  notch,  (b)  IL-2  surface  concentration  in  new  flow  cell. 
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2.1.14  Summary  of  Accomplishments:  Model  Development 

During  the  course  of  this  project,  several  simulation  models  have  been  developed  and 
implemented  in  CFD-ACE+.  Using  these  models,  several  physical  phenomena  have  been 
studied  in  detailed.  These  studies  directly  benefited  the  project  members,  as  well  as  academic 
and  industry  community  at  large.  These  accomplishments  are  discussed  in  detail  next: 

Models  Developed 

Ml.  Particle  dielectrophoreis  models  that  include  conventional  DEP,  traveling  wave  DEP, 
and  electrorotation  [29] 

M2.  AC  electrohydrodynamic  phenomenon  of  electrothermal  flows  [16] 

M3.  DC  electroosmosis  in  systems  with  thick  double  layer  [P23] 

M4.  Slip  boundary  models  for  flow  in  hydrophobic  microchannel  systems  [33,  34] 

M5.  Models  for  EOF  mobility  variations  as  a  function  of  pH  for  fused  silica  substrates [3  5] 

M6.  Unit-cell  based  models  for  highly  concentrated  buffer  solutions.  Models  predict 

following  physical  properties[P23]: 

a.  Electrophoretic  mobility 

b.  Diffusion  coefficient 

c.  Electrical  conductivity 

M7.  Integration  of  biochemistry  models  with  particle  DEP  models  [P45] 

M8.  Joule  heating[36] 

M9.  Super  EU  solver  for  highly  non-net  neutral  system[37] 

MIO.  Generalized  surface  chemistry  models  (beads  and  surfaces)[38] 

Mil.  Biochemistry  database  integrated  with  electrochemistry 

Ml 2.  Hydrogel  models  for  surface  biochemistry[30] 

Ml  3.  Least  square-based  engine  for  extraction  of  kinetic  coefficients [3 8] 

Ml 4.  Rapid  ANN  based  kinetic  coefficient  extraction  algorithm  (this  has  not  been  published 

yet) 

Ml 5.  Python  based  scripts  to  automate  geometric  generation  for  proteomic  chips  (these 

scripts  are  available  from  ESI-CFD,  NA,  Huntsville,  AL) 

The  above  models  been  validated  with  following  analytical  solutions: 

1 .  Joule  heating  in  parallel  plate  [  1 6] 

2.  Particle  DEP  in  periodic  electrode  array[P28] 

3.  Particle  DEP  in  wedge  shaped  (conical)  electrode  system[P28] 

4.  Electrothermal  flow  in  periodic  electrode  array[P25] 
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5.  Ion  transport  near  charged  surface 

6.  EOF  flow  in  arbitrary  zeta  potential  and  double  layer  thickness 

7.  All  surface  chemistry  models  were  validated  against  the  solutions  in  Lu,  et  al.[39] 


After  validation,  the  models  have  been  used  extensively  to  study  various  physical  phenomena, 
and  are  applied  in  the  design  and  development  of  devices.  Table  2.5  summarizes  these  efforts. 


Table  2.5:  Physical  phenomena,  their  models  and  examples  of  the  application  of  the  models. 


Physical  Phenomenon 

Model  Used 

Device  Designed 

1 

Sample  dispersion  under 
non-uniform  zeta  potential 

M3,M5 

2 

Electrothermal  flow 

M2 

Tunable  laser  cavity  sensor  (Fig  1.37) 

3 

Field  Amplified  Sample 
Stacking 

M5,M6 

Sample  stacking  device  (Fig.  1.46,  1.52, 
1.57,  2.4) 

4 

Slip  at  nanoscale 

M4 

5 

Conventional  and  traveling 
wave  DEP 

Ml 

Flow  field  fractionation  devices 

6 

Nanoelectrokinetic 

phenomena 

M3,M6 

Nanofluidic  devices 

7 

Transport  of  analytes  in 
highly  concentrated  buffers 

M6 

8 

Bead  based  immunoassay 

M10,M11 

Notional  immunoassay  based  sensors 

9 

Mass  transport  limitation  in 
biosensor  devices 

M10,M11, 

M12,M15 

Biacore  biosensor 

10 

Effect  of  hydrogel  layer 
(porous  media  flow)  on 
proteomic  chip  performance 

M10,M11, 

M12,M15 

HTS  proteomic  chip 

The  models  enabled  rapid  development  and  demonstration  of  the  devices  mentioned.  In 
particular,  the  models  were  invaluable  in  screening  design  concepts  (HTS  and  tunable  laser 
cavity  sensor)  and  in  data  analysis  (Biacore  biosensor).  In  the  absence  of  the  models,  a  similar 
end-point  may,  if  at  all,  have  been  reached  only  after  expensive  and  time-consuming  iterations. 


Design  guidelines  have  been  generated  for  engineers  and  scientists  who  work  in  this  field.  These 
are  disseminated  to  the  community  via  publication[40],  workshops,  conference  presentation  and 
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tutorials.  Besides,  the  simulation  models  have  been  implemented  in  CFD-ACE+  Version  2002, 
2003  and  2004,  and  have  been  successfully  marketed. 


Commercial  Release  of  Developed  Software  Capabilities:  Both  the  electrokinetic  models  and 
the  biochemical  assay  models  have  also  been  fully  integrated  into  the  commercial  CFD-ACE+ 
code,  along  with  the  development  of  associated  GUI,  user  manuals  and  tutorial  test  cases. 


Master  Classes  for  Advanced  CFD-ACE+  Users:  Sample  test  cases  have  also  been  included 
based  on  these  capabilities  in  the  Master  Classes  (for  advanced  CFD-ACE+  users)  for  the 
biochemistry  and  the  electrokinetic  modules.  These  classes  were  taught  at  the  Annual  CFD- 
ACE+  Users  Conference. 


Multi-physics  Courses:  This  material  (test  cases,  tutorial  examples,  etc.)  is  also  be  included  in 
the  following  two  Biotechnology  industry-specific  multiphysics  courses  being  developed  at 
CFDRC,  which  were  offered  starting  November,  2001  in  various  locations  around  the  country. 
The  two  courses  offered  dealt  with:  (a)  design  of  biochemical  assays  in  microfluidic  systems, 
and  (b)  electrokinetic  and  electrochemistry  applications  in  microfluidic  systems.  The  objective 
of  these  courses  was  to  facilitate  rapid  learning  and  confidence  building  in  existing  and  new 
users  of  the  software  described  here  by  discussing  theoretical  background  (adequate  physical 
models,  data  sources,  and  boundary  conditions),  optimal  computational  strategies,  model  set-up 
protocols  of  CFD-ACE-I-,  and  expected  results  (accuracy,  speed,  uncertainty). 


2.2  Transport  Enhancement:  Electro-thermal  Flow  Modeling  (UCSB) 

AC  electrokinetics  was  used  to  enhance  microfluidic  immuno-sensors,  for  example, 
immunoassays,  in  which  a  ligand  immobilized  on  a  microchannel  wall  specifically  binds 
analyte  flowing  through  the  channel.  These  sensors  can  be  limited  in  both  response  time 
and  sensitivity  by  the  diffusion  of  analyte  to  the  sensing  surface.  The  sensitivity  and 
response  of  these  heterogeneous  immunoassays  may  be  improved  by  using  AC 
electrokinetically-driven  microscale  fluid  motion  to  enhance  antigen  motion  towards 
immobilized  ligands.  Specifically,  the  electrothermal  effect  is  used  to  micro-stir  analyte 
near  the  binding  surface.  Numerical  simulations  of  antigen  in  a  microchannel  flow 
subjected  to  the  electrothermal  effect  show  that  6  Vrms  applied  to  electrodes  near  a  binding 
region  can  increase  binding  in  the  first  few  minutes  by  a  factor  of  seven.  The  effectiveness 
of  electrothermal  stirring  is  a  strong  function  of  the  Damkohler  number.  The  greatest 
binding  enhancement  is  possible  for  high  Damkohler  number,  where  the  reaction  is  limited 
by  diffusion.  Based  on  these  results,  this  technique  is  useful  for  a  large  variety  of 
microfluidic  sensors. 
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2.2.1  Background 


Immunoassays,  which  rely  on  specific  antigen-antibody  binding  for  identification  of  proteins  in  a 
sample,  have  applications  in  both  clinical  laboratories  for  medical  diagnostics  and  treatment 
monitoring,  and  in  research  laboratories  for  highly  multiplexed  testing,  such  as  for  biomarker 
identification.  In  these  cases,  throughput  is  a  key  consideration.  One  factor  that  can  limit  test 
duration  is  diffusion  rate  of  analyte  to  the  reporter.  This  is  particularly  true  for  high  sensitivity 
ELISA  tests  [41].  An  incubation  step  of  minutes  to  hours  allows  diffusion-limited  binding  to 
reach  detectable  levels.  Tests  are  often  performed  at  centralized  labs  where,  despite  a  long  test 
duration,  high  throughput  is  achieved  through  robotics  and  highly  parallel  assays.  However,  if 
the  assay  could  be  moved  from  the  centralized  lab  to  the  point  of  care,  real  time  immunoassay- 
based  diagnostics  could  be  performed.  In  order  for  this  to  be  achieved,  the  test  must  be  faster,  as 
well  as  more  portable,  while  maintaining  high  sensitivity. 

In  response  to  the  needs  for  increasing  throughput,  portability,  and  sensitivity,  new  formats  for 
miniaturized  immunoassays  have  developed  dramatically  in  recent  years.  These  include  spotted 
microarrays,  common  for  “gene  chips”,  and  now  used  for  “protein  chips”  (e.g.  the  ProtoArray 
from  Invitrogen,  Carlsbad,  CA),  and  various  forms  of  Lab-on-a-chip  devices  which  can  perform 
fluid  processing  and  detection  steps  on  a  single  chip  [42].  Small  length  scales  permit  small 
sample  sizes  and  shorter  test  times;  on-chip  sample  preparation  reduces  fluid  handling  steps. 
Though  greatly  aided  by  their  small  length  scales,  these  assays  can  still  be  limited  in  response  by 
diffusion  of  the  analyte  to  an  immobilized  ligand.  In  this  report,  a  method  where  AC 
electrokinetics  is  used  to  enhance  the  performance  of  heterogeneous  assays  is  presented.  Here 
electrothermal  forces  are  used  to  micro-stir  the  analyte  near  a  functionalized  surface,  increasing 
the  rate  of  transport  to  the  surface.  This  addresses  the  need  for  faster  assays  by  offering  a  tool 
that  is  adaptable  to  a  wide  variety  of  assay  configurations  and  reduces  the  incubation  time 
required  to  perform  a  specific  test  while  maintaining  its  sensitivity. 

2.2.2  AC Electrokinetic  Phenomena 

AC  electrokinetics  refers  to  induced  particle  or  fluid  motion  resulting  from  externally  applied 
AC  electric  fields.  DC  electrokinetics  has  been  widely  used  for  lab-on-a-chip  applications  such 
as  electroosmotic  pumping  [43]  and  capillary  gel  electrophoresis  for  DNA  fractionation 
[44], [45].  In  contrast,  AC  electrokinetics  has  received  significantly  less  attention.  AC 
electrokinetics  have  the  advantages  of  (1)  largely  avoiding  electrolysis,  and  (2)  operating  at 
lower  voltages  (1~20  V),  which  is  important  for  portable  systems.  AC  electrokinetics  can  be 
classified  into  three  broad  areas:  dielectrophoresis  (DEP),  electrothermal  forces  and  AC  electro¬ 
osmosis  [16]. 

Dielectrophoresis  takes  advantage  of  a  force  arising  from  differences  in  polarizability  between 
the  particle  and  the  fluid  medium  in  the  presence  of  a  non-uniform  electric  field.  DEP  has  been 
used  to  separate  blood  cells  and  to  capture  DNA  molecules  [46];  [47]  provides  an  overview. 
However,  since  the  DEP  force  scales  with  the  cube  of  particle  radius,  its  effectiveness  for 
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manipulation  of  submicron  molecules  (such  as  10  nm-scale  antigen)  is  limited  in  devices  with 
much  larger  length  scales. 

AC  Electroosmosis  arises  when  the  tangential  component  of  the  electric  field  interacts  with  a 
field-induced  double  layer  along  a  surface.  A  bioprocessor  has  been  developed  at  UCLA  [48] 
for  the  concentration  of  bioparticles  including  bacteria  and  ^-phage  DNA.  This  device  relies  on 
the  balance  between  electroosmotic  flow  and  DEP  force  on  suspended  particles.  In  another 
application,  Bazant  and  Squires  [49]  develop  AC  electroosmosis  for  pumping  by  inducing  a 
double  layer  around  metal  posts  that  are  asymmetric  either  in  shape  or  in  surface  properties.  AC 
Electroosmosis  is  only  effective  when  there  is  a  substantially  deep  induced  double  layer,  that  is, 
for  low  conductivity  solutions  and  low  applied  field  frequencies.  For  example,  in  an  aqueous 
saline  solution  with  an  electrical  conductivity  of  a  =  2  x  10'^  S/m,  it  is  predicted  that  AC 
electroosmosis  is  not  important  for  frequencies  />  100  kHz  [50]. 

Transport  enhancement  of  small  proteins  may  be  most  successful  through  electrothermally 
driven  flow  (ETF).  A  non-uniform  electric  field  produces  uneven  Joule  heating  of  the  fluid, 
which  gives  rise  to  nonuniformities  in  conductivity  and  permittivity.  These  interact  with  the 
electric  field  to  generate  flow,  often  in  circulating  patterns.  Unintentionally-induced 
electrothermal  flow  in  a  traveling  wave  DEP  cell-sorter  was  studied  numerically  [51],  and  found 
to  be  quite  sensitive  to  thermal  boundary  conditions  such  as  channel  material,  and  applied 
temperature  differential  boundary  conditions.  Characteristic  swirling  flow  patterns  can  be  used 
to  circulate  suspended  molecules  past  the  binding  region,  thereby  providing  more  binding 
opportunity  for  the  suspended  molecules. 

2.2.3  Experimental  Micro-PIV  Measurements 

With  the  objective  of  exploring  two-dimensional  fluid  and  particle  motion  in  a  non-uniform 
electric  field,  a  cavity  was  fabricated  with  sidewall  electrodes  (Figure  2.37).  The  device  was 
sandwiched  between  glass  slides  and  filled  with  a  solution  of  0.05  M  KCl  and  fluorescent 
polystyrene  tracer  particles.  When  the  electrodes  were  driven  at  E  =  7  Vrms  and /=  200  kHz,  a 
fluid  circulation  pattern  was  observed.  The  two-color  micron-resolution  particle  image 
velocimetry  (p-PIV)  technique  presented  in  Wang  e(  al.  [52]  was  used  to  measure  fluid  motion. 
These  experiments  are  described  in  detail  in  [53]  and  are  summarized  here  to  provide  motivation 
for  the  numerical  simulations.  The  measured  fluid  velocity  field  is  shown  in  Figure  2.38.  The 
applied  field  induces  a  circulating  flow,  with  a  characteristic  velocity  of  about  100  pm/s.  When 
the  driving  voltage  is  varied  between  3-7  Vrms,  the  flow  pattern  remains  self  similar,  but  the 
velocity  magnitude  varies  with  voltage  to  the  4*  power,  which  is  characteristic  of  electrothermal 
flow  [16]  (data  not  shown).  The  treatment  of  electrodes  as  isothermal  is  appropriate  for 
electrodes  of  sufficient  thickness  relative  to  length. 

Gradients  in  temperature  produce  gradients  in  electrical  permittivity  and  conductivity  in  the 
fluid.  For  water,  the  dependence  of  electrical  conductivity  on  temperature  is  (7/cr)  (dcr/ffT)  = 
+2%  and  the  dependence  of  electrical  permittivity  on  temperature  is  {1/s)  {ds/dT)  =  -0.4%  per 
degree  Kelvin  [54].  These  variations  in  electric  properties  produce  gradients  in  charge  density 
and  perturb  the  electric  field.  Assuming  the  perturbed  electric  field  is  much  smaller  than  the 
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applied  electric  field,  and  that  advection  of  electric  charge  is  small  compared  to  conduction,  the 
time-averaged  electrothermal  force  per  unit  volume  for  a  non-dispersive  fluid  can  be  written  as 
[16]: 
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where  r  =  sja-  is  the  charge  relaxation  time  of  the  fluid  medium  and  the  incremental  temperature- 
dependent  changes  are: 
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Figure  2.37:  Test  device  used  to  measure  2D 
flow  created  by  an  inhomogeneous  AC  electric 
field. 

The  device  is  sandwiched  between  glass,  and  the 
cavity  is  filled  with  0.05  M  KCl  solution  seeded 
with  fluorescent  polystyrene  particles.  An  AC 
driving  signal  produces  an  electric  field  in  the 
fluid  which  generates  fluid  motion. 


Figure  2.38:  Fluid  velocity  field  measured  in  the 
device  shown  in  Figure  2.37.  Images  of  the 
fluorescent  particles  moving  in  the  fluid  are 
captured  at  discrete  times.  Micro-Particle  Image 
Velocimetry  is  used  to  determine  particle  velocity 
field.  Data  from  two  different  size  particle 
velocity  fields  are  combined  to  estimate  the  fluid 
velocity.  The  electrodes  are  driven  at  7  Vrms  and 
200  kHz.  The  AC  field  produces  two  counter¬ 
rotating  vortices  of  characteristic  velocity  100 
pm/s. 
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The  first  term  on  the  right  hand  side  of  Equation  (2.43)  is  the  Coulomb  force,  and  is  dominant  at 
low  frequencies,  «=  27if  The  second  term  is  the  dielectric  force,  and  is  dominant  at  high 
frequencies.  The  crossover  frequency  scales  inversely  with  the  charge  relaxation  time  of  the 
fluid  [16].  For  example,  an  aqueous  solution  with  conductivity  (j=  10'^  S/m  has  a  crossover 
frequency  around /=  10  MHz. 

The  electrothermal  force  shown  in  Equation  (2.43)  is  a  body  force  on  the  fluid.  The  motion  of 
the  fluid  can  be  determined  by  solving  the  Stokes’  equation  for  zero  Reynolds  number  fluid  flow, 
such  that: 


0  =  -  Vp  +  pV  ^  u  +  F^j. 

r 

V  gu  =  0 

(2.45) 

1 

where  u  is  the  fluid  velocity,  p  is  the  pressure  in  the  fluid,  and  p  is  the  dynamic  viscosity  of  the 
fluid.  Figure  2.39  shows  the  resulting  velocity  field.  With  an  applied  voltage  of  F  =  5  Vrms,  the 
velocity  of  the  electrothermal  flow  is  on  the  order  of  100  pm/s  and  is  characterized  by  a  pair  of 
counter  rotating  vortices.  This  velocity  field,  simulated  for  V  =  5  Vrms,  closely  matches 
experimentally  measured  velocity  for  V  =  1  Vrms  (Figure  2.38).  Out  of  plane  heat  transfer, 
unique  to  the  thin  experimental  device,  is  not  accounted  for  in  the  2D  model,  and  may  cause  this 
difference.  The  model  was  run  with  several  grid  resolutions,  and  determined  to  be  grid- 
independent. 

The  total  power  consumption  was  calculated  by  integrating  the  Joule  heating  term,  cE  ,  over  the 
flow  domain.  For  this  particular  model  (ow  =  0.66  S/m;  applied  potential  V  =  5  Vrms  at /=  200 
kHz),  and  for  an  out  of  plane  cavity  depth  d  =  200  pm,  the  predicted  power  dissipated  is  4.3 
mW. 
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Figure  2.39;  Simulation  of  electrothermally-driven  fluid  motion  with  an  applied  voltage  of  5  Vrms. 

The  velocity  of  the  electrothermally-driven  flow  is  of  order  100  pm/s  and  is  characterized  by  a  pair 
of  counter  rotating  vortices.  This  flow  pattern  matches  the  velocity  field  that  was  experimentally 
measured  using  micro-PIV  (see  Figure  2.38)  with  an  applied  voltage  of  7  Vrms- 
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Figure  2.40:  Numerical  simulation  of  temperature  distribution  created  by  Joule  heating. 

The  maximum  temperature  rise  of  2.9  K.  The  sidewall  electrodes  are  separated  by  25  pm,  and  an 
applied  voltage  of  5  Vrms  are  applied  to  the  electrodes. 


56 


500  X  600  ^im 


(a)  Diffusion  only 


depleted  concentration 
surrounds 

functionalized  surface 


50  pm  ^ 

functionalized  surface 


(b)6V, 


depleted  concentration  better 
mixed  throughout  domain 


z 


electrodes  with 
25  pm  gap 


T 


functionalized 
surface  on  gold 
electrode 


CC„-’ 

11.0 

0.9 

0.8 

I 

0.6 

0.5 

0.4 


0.3 


0.2 

0.1 

0 


Figure  2.41 :  Simulation  of  concentration  in  a  microcavity  after  100  sec,  subject  to  binding  on  the 

functionalized  surface  at  the  wall. 


In  (a),  where  mass  transport  is  through  diffusion  only,  a  depleted  cloud  (black)  surrounds  the 
functionalized  surface,  indicating  a  depletion  of  analyte.  In  (b),  electrothermal  flow  circulates 
the  depleted  cloud  throughout  the  cavity,  and  circulates  fresh  analyte  past  the  functionalized 
surface,  yielding  a  higher  effective  concentration  close  to  the  functionalized  surface. 

2.2.4  Effect  on  Binding:  Simulation  of  Microcavity 

The  highest  flow  velocities  are  predicted  directly  over  the  inside  edges  of  the  electrodes. 
Because  of  the  high  transport  rate  here,  one  expects  this  to  be  an  optimum  sensor  location.  The 
gold  electrode  surface  provides  a  convenient  surface  for  antibody  immobilization  through  thiol 
linkers.  Therefore,  the  effect  of  this  flow  pattern  on  the  binding  response  of  an  assay  in  which 
antibody  has  been  immobilized  along  a  short  length  of  the  electrode  near  the  electrode  gap  is 
investigated  (Figure  2.41).  For  these  microstirring  enhanced  binding  simulations,  a  slightly 
higher  potential  is  applied,  6  Vrms»  than  in  the  previous  simulations  (5  Vrms),  for  faster  stirring. 
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The  convective  scalar  equation  describes  the  suspended  concentration  c(x,y,t)  of  antigen  within 
the  microchannel: 

—  +  i-Vc  =  DV^c 

dt 

(2.46) 


where  u  is  the  fluid  velocity  and  D  the  diffusivity  of  the  antigen.  An  initial  concentration  in  the 
cavity  Co  is  depleted  through  binding  at  the  wall.  Following  the  model  given  by  Myszka  [55],  the 
rate  of  binding  at  the  wall  for  a  first  order  reaction  is  k^^  c^(Rj.-B),  where  is  the  on  rate 

constant,  Rr  is  the  receptor  concentration,  B  is  the  bound  antigen  concentration,  and  Cw(x)  is  the 
suspended  concentration  of  antigen  along  the  wall.  The  rate  of  dissociation  is  kof/B  ,  where  kgff 
is  the  off  rate  constant.  The  time  rate  of  change  of  antigen  bound  to  the  immobilized  antibodies, 
dBjdt ,  is  equal  to  the  rate  of  association  minus  the  rate  of  dissociation: 


—  =  k  c  (R^  —  B)  —  k  ^  B 

on  w  \  T  /  off  ' 

(2.47) 

The  rate  of  antigen  binding  to  immobilized  antibodies  must  be  balanced  by  the  diffusive  flux  of 
antigen  at  the  binding  surface,  y  =  0,  such  that: 


(2.48) 


dB 

dt 
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Figure  2.42  :  Normalized  binding  curves  for  non-enhanced  (0  Vnns)  and  electrothermally  enhanced  (6 

Vrms,  12  Vnns)  transport. 


B  is  the  amount  of  analyte  bound,  which  is  normalized  by  the  concentration  of  immobilized 
antibodies,  Rj .  The  differences  in  the  curves  show  an  increase  in  binding  rate  which  yields  a 
factor  of  7  higher  binding  for  an  applied  voltage  of  6  Vrms  and  a  factor  of  14  higher  binding  after 

100  seconds  for  12  Vrms- 


Equations  (2.46),  (2.47)  and  (2.48)  are  solved  with  an  immobilized  antibody  concentration  Rj  = 
3.3  X  lO  "  M  m  (i.e.  2  x  lO'^  molecules/m^)  and  an  initial  suspended  antigen  concentration  of 
10''*^  M.  Under  the  standard,  passive  format  where  there  is  no  micro-stirring  and  binding  locally 
depletes  the  suspended  antigen  concentration.  Because  diffusion  is  the  only  transport 
mechanism,  the  depleted  region  surrounding  the  sensor  surface  grows  with  time  and  reduces  the 
rate  of  binding.  When  the  functionalized  surface  is  located  near  the  center  of  an  electrode  pair 
(see  Figure  2.44),  and  the  electrodes  are  energized  at  F  =  6  Vrms  and  /  =  200  kHz, 
electrothermally-driven  circulation  redistributes  the  depleted  concentration  throughout  the 
domain,  and  the  sensor  effectively  sees  higher  suspended  analyte  concentration,  resulting  in  a 
higher  binding  rate.  Figure  2.43  shows  that  an  applied  voltage  of  F  =  6  Vrms  produces  a  7-fold 
increase  in  bound  antigen  compared  to  the  0  Vrms  (passive)  case.  For  an  applied  voltage  of  F  = 
12  Vrms,  the  increase  jumps  to  14-fold.  The  associated  velocities  (wmax  ~  0.4  -  6  mm/s)  vary 
with  voltage  to  the  4*  power,  which  is  characteristic  of  electrothermal  flow.  These  velocities  are 
characteristic  of  other  types  of  AC  electrokinetic  flows.  For  example,  Bazant  and  Squires  predict 
a  0.7  mm/s  velocity  generated  by  charge-induced  AC  electroosmosis  [49].  To  be  conservative, 
this  work  will  focus  on  the  6  Vrms  case,  for  which  the  numerical  simulations  suggest  a  factor  of 
7.2  increase  in  binding  rate.  Higher  binding  rates  may  be  achieved  with  higher  voltages.  The 
following  parameters  were  used  in  the  numerical  simulation:  solution  conductivity  ow  =  066  S 
m'';  on  and  off  rate  constants  kon  =  10^  M''s''  and  kof^W^  s'';  immobilized  receptor 
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Figure  2.43  :  Velocity  field  in  microchannel  showing  effect  of  electrothermal  force  on  channel  flow. 

Inlet  mean  velocity:  67  (j,m/s.  Electrothermal  forces  from  a  6  Vrms  electric  field  generate 
velocities  which  are  on  the  order  of  300  (rm/s. 


concentration  Rr  =3.3  x  10‘”  Mm;  diffusivity  D  =  10'"  m^s"';  initial  suspended  antigen 
concentration  Co  =  0. 1  nM.  The  elements  in  the  2-D  model  are  triangular  with  sizes  ranging  from 
1  pm  near  the  electrode  gap  to  20  pm  far  from  the  gap.  At  this  mesh  size,  the  solution  was 
shown  to  be  mesh-independent.  The  boundary  condition  at  the  functionalized  surface  (Eq.  2.48) 
is  solved  in  a  separate  one  dimensional  domain  and  coupled  to  the  2D  model.  In  order  to 
facilitate  the  solution  convergence,  which  can  be  difficult  for  low  diffusivities,  artificial 
streamline  diffusion  is  added  to  the  model.  It  was  verified  that  this  does  not  affect  the  results: 

the  average  bound  concentration  B(t  =  100^)  agrees  within  0.3%  for  cases  with  and  without 
streamline  diffusion  (data  not  shown). 

2.2.5  Flow-through  Simulations 

The  team  investigated  the  binding  enhancement  from  micro-stirring  in  microchannels  with 
pressure-driven  flow.  In  these  sensors,  the  pressure-driven  flow  continually  exposes  the  sensor 
to  fresh  analyte,  and  so  one  does  not  expect  that  the  binding  enhancement  through  micro-stirring 
will  be  as  significant  as  for  non-pressure  driven  flow.  The  microchannel  is  also  smaller  in  these 
flow-through  simulations  to  better  represent  microscale  flow-through  channels  (e.g.  MEMS 
based)  which  often  have  smaller  fluid  volume  than  static  immunoassays  (eg.  microarray).  The 
channels  are  100  pm  in  height,  so  the  Damkohler  number  is  reduced  by  a  factor  of  5,  as 
compared  with  the  same  assay  in  the  larger  microcavity  simulations.  For  these  simulations,  a 
parabolic  velocity  inlet  boundary  condition  is  prescribed.  Electrothermal  forces  distort  the 
parabolic  flow  around  the  electrode  gap  to  produce  a  circulating  flow  as  shown  in  Figure  2.43. 
For  this  case  (average  inlet  velocity  u  =  61  pm/s;  V  =  6  Vrms)  the  electrothermally  generated 
velocity  (~300  pm/s)  is  large  compared  with  the  average  base  flow  (67  pm/s).  In  the  convection- 
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diffusion  model,  zero  initial  concentration  is  prescribed  throughout  the  microchannel.  A 
concentration  of  Co  =  0. 1  nM  is  introduced  at  the  inlet  (left  hand  side)  of  the  channel  for  t  >  0. 
For  the  passive  case,  where  there  is  no  applied  voltage,  the  base  flow  is  parabolic,  and  analyte  is 
transported  downstream  most  rapidly  at  the  channel  center.  This  is  shown  in  Figure  2.45a,  where 
the  highest  concentrations  of  analyte  extend  through  the  center  of  the  channel.  When  an  electric 
potential  is  applied,  the  electrothermally-induced  motion  redirects  the  concentration  profde 
(Figure  2.45b).  The  optimum  position  of  the  functionalized  surface  is  near  the  leading  edge  of 
the  downstream  electrode. 


In  this  case,  where  convective  transport  from  pressure-driven  flow  through  the  microchannel  is 
important,  the  Peclet  number  must  be  considered  as  well  as  the  Damkohler  number.  The  Peclet 
number  is  defined  as  Pe  =  uhD ,  which  is  the  ratio  of  convective  to  diffusive  transport. 
Figure  2.45  shows  the  binding  enhancement  factor  Be  for  t  =  100  s,  for  various  Peclet  numbers, 
and  as  a  function  of  Da.  As  in  the  non-flow-through  format  ( 

Figure  2.46),  the  binding  enhancement  factor  Be  has  a  strong  dependence  on  Da.  In  this  flow¬ 
through  case,  however,  there  is  an  additional  strong  dependence  on  Pe. 

Figure  2.46  indicates  that  the  binding  enhancement  factor  Be  is  lower  for  larger  Pe,  that  is,  for 
faster  flows.  While  electrothermal  micro-stirring  may  not  be  able  to  increase  reaction  rates  in 
channels  with  high  flow  rates,  it  may  provide  more  efficient  capture  of  antigen,  thereby  allowing 
reduced  flow  rates.  In  cases  where  a  high  flow  rate  has  been  required  to  achieve  sufficiently  high 
binding  rates  in  a  transport-limited  regime,  the  addition  of  electrothermal  stirring  to  the  channel 
allows  the  flow  rate  to  be  reduced  while  maintaining  the  same  high  binding  rate.  This  effectively 
reduces  required  sample  volume. 


(a)0V 
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Figure  2.44:  Concentration  plots  of  electrothermally  modified  channel  flow  at  t  =  5  sec  with  applied 

voltages  of  (a)  0  Vrms  and  (b)  6  Vn„s. 


With  optimal  size  and  placement  of  electrodes,  the  electrothermal  eddies  can  be  engineered  to 
span  the  width  of  the  channel,  as  is  the  case  here,  for  a  100  micron  channel.  High  flow  rates  and 
concentration  gradients  and  therefore  an  increase  in  diffusive  flux  in  the  vertical  direction 
directly  over  the  inside  edges  of  the  downstream  electrode  indicate  a  favorable  location  of  the 
functionalized  surface. 
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Figure  2.45;  Simulation  results  show  binding  enhancement  factor,  Be,  increases  with  increasing 
Damkohler  number  Da  as  in  the  static  case,  but  decreases  with  increasing  Peclet  number,  that  is,  faster 

pressure-driven  flow. 


2.2. 6  Binding  Enhancement  Summary 

Electrothermally-generated  fluid  motion  was  investigated  experimentally  and  numerically. 
Experimental  measurements  of  AC  electrokinetically-driven  velocity  fields  suggest  that 
electrothermal  forces  are  more  important  than  AC  electroosmosis  for  generation  of  fluid  motion 
in  high  conductivity  (a  =  0.66  S/m)  buffers,  which  are  applicable  for  many  bioassays. 
Numerical  simulations  were  conducted  to  predict  electrothermally-generated  fluid  motion.  The 
model  predicts,  for  5  Vrms  applied,  a  temperature  rise  on  the  order  of  2.9  K,  which  results  in 
velocities  on  the  order  of  100  pm/s.  The  numerical  model  flow  pattern  closely  matches 
experimental  measurements  with  an  applied  voltage  of  F  =  7  Vrms-  The  total  power  consumption 
of  the  electrothermal  flow  was  4.3  mW  for  a  200  pm  deep  cavity. 

This  model  was  used,  together  with  a  first  order  reaction  model,  to  predict  binding  enhancement 
in  heterogeneous  reactions,  for  example,  immunoassays.  The  model  suggests  that 
electrothermally-generated  micro-stirring  can  increase  the  binding  rate  for  non-flow-through 
assays,  such  as  a  microarray  or  microtiter  plate,  by  a  factor  of  seven,  for  6  Vrms  applied  potential. 
A  series  of  simulations  with  different  parameters  show  that  the  binding  enhancement  depends 
strongly  on  the  Damkohler  number  (Da).  The  binding  rate  is  increased  by  a  factor  of  7.3  with 
electrothermal  micro-stirring,  for  Da  =  lO"^.  However,  the  binding  rate  is  only  increased  by  a 
factor  of  3  for  Da  ~  100,  and  little  increase  in  binding  rate  is  observed  for  Da  ~  1 . 

Electrothermally-generated  micro-stirring  can  also  improve  binding  rates  for  flow-through 
assays,  for  example  in  BioMEMS  or  lab-on-a-chip  devices.  In  these  cases,  the  binding 
enhancement  factor  depends  on  the  Peclet  number  as  well  as  the  Damkohler  number.  Higher 
Peclet  numbers,  which  correspond  to  higher  flow  rates,  yield  lower  binding  enhancement 
through  electrothermal  micro-stirring.  For  example,  with  a  Damkohler  number  of  3  x  10^,  and  a 
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Peclet  number  of  100,  a  factor  of  nearly  8  improvement  is  possible.  However,  with  a  Peclet 
number  of  1000,  the  expected  improvement  is  reduced  to  a  factor  of  about  3.  The  numerical 
simulations  reported  here  indicate  that  electrothermal  micro- stirring  can  be  a  useful  technique  for 
increasing  binding  rates  in  heterogeneous  assays,  particularly  for  diffusion-limited  reactions. 

2.3  Design  and  Fabrication  of  Electrokinetic  Systems  (ACLARA  BioSciences) 


2.3.1  Plastic  Microfluidic  Device  Fabrication 

The  design  and  fabrication  of  microfluidic  devices  in  plastic  has  tremendous  potential  impact  in 
the  commercial  endeavor  of  pharmaceutical  drug  discovery,  where  such  assays  as  multiplexed 
gene  expression  analysis,  multiplexed  proteomics  experiments,  and  multiplexed  enzyme  assays 
have  been  pioneered  by  ACLARA  and  other  fluidics  companies  and  research  groups  to  lower  the 
cost,  enhance  the  performance,  and  speed  the  process  of  these  bioassays.  To  optimize  such 
measurements,  a  wide  range  of  designs  and  microfluidic  configurations  have  been  designed  and 
prototyped.  Plastic  microfluidics  also  are  beginning  to  have  significant  impact  in  clinical 
diagnostics,  where  disposability  is  a  must  and  functional  integration  to  save  size,  cost,  and 
sample  volume  while  improving  reliability  are  paramount.  The  Stanford  DARPA  SIMBIOSYS 
Program  pushed  the  frontiers  of  such  parameters  as  microchannel  aspect  ratio,  thus  helping 
ACLARA  to  define  the  limits  of  low-cost  plastic-to-plastic  card  technologies.  The  multiplexed 
assays  mentioned  above  were  also  used  to  identify  biological  pathogens  under  DARPA 
BioFLIPS  Program,  a  synergistic  effort  with  the  SIMBIOSYS  Program. 

In  addition,  ACLARA  developed  expertise  and  demonstrated  technology  in  several  key  areas 
that  contributed  to  rapid  progress  in  the  program: 

•  low-defect  glass  mastering  of  fluidic  chips 

•  RIE-based  mastering  (using  the  Stanford  facility,  and  commercial  vendors)  for  variable 
aspect  ratios  and  unique  geometries 

•  prototyping  of  moderate  volumes  (5  -  50)  of  new  designs 

•  thermal  lamination  of  plastic-to-plastic  chips  to  maintain  channel  definition 

•  pressure-sensitive  adhesive  channel  sealing  to  enable  in-chip  reagent  storage 

•  multilevel  channel  patterns  for  complex  fluidic  architectures 

•  integrated  heaters  for  on-chip  temperature  control/monitoring 

•  integrated  electrodes  for  on-chip  drive  of  electrokinetic  processes 

•  laser-induced  fluorescence  detection  in  plastic  chips 

•  high-resolution  CE  separations  in  plastic  microchannels 

•  highly  multiplexed  bioassays 
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2.3.2  Activities 


Microfabrication  of  Masters  with  Complex  Geometries  using  Silicon  RIE,  SU-8  Patterning,  and 
Glass  Etching. 

ACLARA  developed  a  new  method  of  ITP  (isotachophoresis,  or  “stacking”)  for  which  simulated 
results  have  been  reported  [56],  This  new  ITP  methodology  offers  high  concentration  ratios  and 
simpler  electrical  control  relative  to  conventional  ITP.  In  general,  these  ITP  designs  include 
geometric  features  requiring  higher  aspect  ratios,  different  angles,  and  better  pattern  resolution 
than  attainable  with  wet  etching  of  glass  plates  (ACLARA’s  best-developed  mastering 
technology).  Further  development  of  the  RIE  process  was  completed,  and  optimal  etching 
parameters  were  selected  to  provide  non-reentrant  sidewall  profdes  with  minimal  sidewall 
roughness  and  micromasking  defects.  A  total  of  six  silicon  masters  were  etched  with  the  point 
concentrator  design  and  evaluated.  An  SEM  photo  of  a  Si  master  and  optical  micrograph  of  a 
molded  plastic  part  from  such  a  master  are  shown  in 
Figure  2.46. 


Figure  2.46:  Left:  SEM  image  showing  reactive  ion  etched  silicon  master  for  point-concentrator  device 

(close-up  of  a  channel  intersection). 

Right:  example  of  a  plastic  part  molded  with  an  electroform  mold  tool  that  was  produced  from  a  reactive 
ion  etched  silicon  master,  sawed  to  show  the  cross-section. 


This  particular  stacking  method,  called  the  “point  concentrator”,  relies  on  a  distortion  of  the 
electric  field  in  the  injection  region  caused  by  a  geometrical  feature  protruding  into  the  channel 
in  the  injection  offset.  The  result  of  the  distortion  to  the  electrical  field  is  an  increased 
concentration  of  ionic  species  in  the  region  of  the  point  concentrator.  The  concentrated  sample 
plug  is  then  injected  and  separated,  to  achieve  a  higher  observed  signal- to-noise  ratio  at  the 
detector.  Point  concentrator  designs  were  fabricated  with  silicon  RIE  and  SU-8  mastering 
processes.  Both  approaches  produced  channel  dimensions  with  close  to  1:1  aspect  ratios. 
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Electroforms  were  successfully  generated  from  these  masters,  and  were  used  to  mold  plastic 
prototype  cards. 


A  set  of  experiments  designed  to  evaluate  the  anticipated  signal  enhancement  relative  to  a 
control  pattern  with  identical  channel  dimensions  and  lengths  was  performed.  Separations  were 
attempted  using  both  1)  40  nM  fluorescein  in  25  mM  HEPES  and  25  mM  NaCl,  and  2)  four 
ACEARA  eTag  reporters[57]  in  lx  buffer  with  12.5  mM  MgCli.  In  both  cases,  a  running  buffer 
of  1%  PEO  in  25  mM  HEPES  was  used,  and  the  separation  was  detected  1  cm  downstream  from 
the  injector.  Injections  were  run  using  both  pinch/pullback  and  floating/stacking  scenarios,  with 
results  shown  in  Figure  2.47. 
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Figure  2.47;  Left:  Separation  performance  of  point  concentrator,  compared  to  control  chip  without  this 

feature,  using  fluorescein. 

Right:  Separation  performance  of  point  concentrator,  compared  to  control  chip  without  this  feature,  using 
4  eTag  reporters  (molecular  tags  of  different  mobilities  that  are  used  in  genomic  and  proteomic 

multiplexed  assays). 


Modeling  results  suggested  that  a  signal  enhancement  of  around  2x  should  be  expected  from  the 
introduction  of  the  point  concentrator  feature  described  above.  However,  for  both  fluorescein 
and  eTag  reporters,  signal  enhancement  was  insignificant  relative  to  the  control  pattern.  The 
likely  cause  of  this  failure  has  been  identified  as  a  rounded  peak  on  the  triangular  point 
concentrator,  which  stems  from  the  difficulty  of  microfabricating  an  arbitrarily  sharp  peak. 
Modeling  results  were  achieved  using  an  idealized  triangular  geometry,  while  the  feature  on  the 
plastic  device  was  slightly  rounded  at  the  feature  peak.  As  shown  in 

Figure  2.48,  both  the  RIE  and  SU-8  processes  resulted  in  slightly  rounded  point  concentrators 
and  channel  comers.  This  is  most  likely  due  to  the  wet  etching  of  the  oxide  mask  in  the  RIE 
case,  and  the  nature  of  photoresist  development  in  the  SU-8  case.  Follow-up  adjustments  to  the 
model  suggested  that  the  effect  of  the  point  concentrator  diminishes  rapidly  as  its  geometry  is 
rounded. 
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To  improve  point  concentrator  performance,  efforts  to  improve  the  sharpness  of  the  point 
concentrator  feature  were  made.  In  particular,  an  isotropic  glass  etching  approach  was 
investigated.  With  isotropic  etching,  it  is  possible  to  obtain  sharp  convex  comers  with  an 
appropriately  designed  etch  mask.  A  variety  of  potential  point  concentrator  features  have  been 
considered  in  addition  to  the  simple  triangular  feature  described  above.  A  transparency 
photomask  was  created  as  the  first  step  toward  quickly  mastering  and  evaluating  a  variety  of 
potential  point  concentrator  geometries.  The  most  promising  features  from  this  mastering  test 
were  incorporated  into  a  glass  photomask.  With  improved  mastering  techniques,  it  is  hoped  that 
both  the  triangular  point  concentrator  and  other  potential  geometries  can  be  demonstrated  to 
result  in  the  signal  increase  suggested  by  the  original  modeling  results. 


Figure  2.49  shows  an  example  of  a  triangular  point  concentration  feature  as  it  appears  in  an 
AutoCAD  design  file,  on  a  transparency  photomask,  and  in  the  etched  glass  master.  This 
triangular  feature  extends  across  the  full  width  of  the  photomask  line.  Using  these  results,  an 
optimized  design  with  appropriate  line  widths  was  used  to  obtain  a  chrome  photomask.  A 
number  of  different  designs  were  considered,  including  triangle-shaped  features  extending  into 
the  channel,  conical  sharp  points  in  the  center  of  the  channel,  and  raised  sharp  ridges  ranning 
perpendicular  or  parallel  to  the  channel.  Figure  2.50  shows  an  example  of  the  sharp  features 
possible  with  isotropic  glass  mastering  techniques. 
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Figure  2.48:  Electroform  (left)  and  molded  plastic  part  (right)  from  a  reactive  ion  etched  silicon  master, 
showing  rounded  channel  corners  and  point  concentrator  features. 


Figure  2.49:  Triangular  point  concentration  feature  across  full  photomask  line  width. 


(a)  CAD  layout;  (b)  transparency  photomask;  (c)  optical  micrograph  of  resultant  etched  glass  master;  (d) 
rotated  perspective  optical  micrograph  of  etched  glass  master. 
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Figure  2.50:  Scanning  electron  micrograph  (SEM)  showing  features  created  with  isotropic  glass  etching 

techniques. 


Fabrication,  inspection,  and  delivery  of  plastic  microfluidic  chips. 

A  number  of  plastic  chips  of  the  DS-1  and  DS-2  chip  designs  were  fabricated,  inspected,  and 
delivered  to  Stanford.  The  DS-2  design  demonstrated  the  production  of  20  pm-deep  channels  in 
thermoplastic,  an  advance  over  typical  plastic  fluidic  cards.  Table  2.6  shows  the  cards  delivered 
to  Stanford. 
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Table  2.6:  Microfluidic  card  deliveries. 


Card  #  (ETD)  \Order 

Design 

Suborder 

1 

Suborder 

2 

Suborder 

3 

Suborder 

4 

Suborder 

5 

DSl  Shallow 

Depth:  20  pm 

Width:  50  and  200  pm 

15 

(03/22/02) 

16 

(05/08/02) 

15 

(07/15/02) 

15 

(08/18/02) 

/ 

DSl  Deep 

Depth:  40  pm 

Width:  90  and  240  pm 

15 

(04/03/02) 

13 

(06/18/02) 

15 

(07/15/02) 

15 

(08/31/02) 

/ 

DS2 

Depth:  20  pm 

Width:  50  pm 

10 

(04/26/02) 

10 

(07/01/02) 

10 

(07/30/02) 

10 

(09/15/02) 

10 

(09/30/02) 

A  second  DS-2  electroform  was  obtained  in  May  2002,  alleviating  the  defect  problems  created 
during  the  de-molding  process  by  using  the  original  DS-2  eform.  Plastic  cards  were  molded 
from  this  new  eform  and  the  image  and  profilometer  measurement  of  the  channel  showed  that 
the  lip  defect  was  eliminated,  resulting  in  improved  channel  replication. 

ACLARA  developed  a  non-destructive  evaluation  technique  for  examining  cross-sectional 
profiles  of  masters  and  electroforms.  The  process  involves  casting  a  layer  of  silicone  rubber 
(polydimethyl  siloxane,  PDMS)  onto  the  substrate  bearing  the  microchannel  features  of  interest; 
curing  the  rubber;  peeling  it  from  the  substrate;  cross-sectioning  with  a  sharp  metal  blade;  then 
examining  the  cross  section  by  optical  or  electron  microscopy.  Figure  2.50  shows  a  channel 
cross  sectional  image  obtained  using  this  method. 


2.3.3  Simulation  and  Demonstration  of  Point  Concentrator  Device  (ACLARA  BioSciences 
and  CFD  Research  Corp.) 

Improvement  of  the  limits  of  detection  in  electrophoretic  separations  can  enable  more  accurate 
and  sensitive  medical  diagnostics,  forensic  measurements,  and  other  bioassays  with  ever  smaller 
sample  quantities.  The  early  diagnosis  of  cancer,  for  example,  can  be  more  effective  as  a 
consequence  of  the  capability  to  detect  fewer  copies  of  a  particular  gene  or  a  smaller 
concentration  of  key  proteins. 

In  the  floating-stacking  approach,  concentrations  of  salts  and  ionic  mobilities  are  selected,  and 
potentials  on  the  various  electrodes  are  controlled  (including  being  left  “floating”),  such  that 
isotachophoresis  (ITP,  or  “stacking”)  conditions  are  established.  Simulations  showed  that  this 
method  can  decrease  sample  loss  and  improve  sensitivity  relative  to  more  standard  methods. 
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The  floating-stacking  experiments  were  conducted  on  ACLARA’s  cTag™  reporter  chemistries  (a 
series  of  functionalized  fluorescein  molecules,  each  with  a  slightly  different  electrophoretic 
mobility;  see  [57])  to  determine  whether  a  large  set  of  multiple  species  behave  in  the  same 
manner  as  the  single  component  simulation.  The  mixture  consisted  of  12  eTag  reporters  in  25 
mM  HEPES  buffer  used  with  a  plastic  chip  with  a  standard  double-T  injector,  3 -mm  offset,  and 
with  the  injection  and  separation  voltage  and  timing  given  in  Table  2.7  and  Table  2.8. 


Table  2.7:  Injection  voltage  setup  and  timing  for  floating/stacking. 


V1 

V2 

V3 

V4 

Time  (sec) 

Inj 

float 

500 

float 

0 

12 

Sep 

0 

float 

4200 

float 

Table  2.8:  “Pinch-and-pullback”  voltage  and  timing  setup. 


V1 

V2 

V3 

V4 

Time  (sec) 

Inj 

0 

500 

0 

0 

12 

Sep 

0 

380 

4200 

380 

In  the  resulting  simulated  electropherograms,  each  of  the  peak  areas  was  integrated  and  the  areas 
for  the  pinch-and-pullback  method  were  compared  to  the  floating/stacking  method.  These  results 
are  summarized  in  Table  2.9,  and  indicate  that  there  is  a  near-uniformity  in  the  ratio  of  areas  for 
floating/stacking  vs.  pinch  &  pull-back  for  each  of  the  12  eTag  reporters.  This  proves  that  the 
electrophoretic  “loading”  of  the  various  eTag  reporters  is  accomplished  in  a  “democratic” 
fashion  (without  bias);  i.e.,  the  ratio  of  eTag  reporter  concentrations  at  the  detector  is  the  same  as 
the  ratio  of  concentrations  in  the  original  sample.  This  critical  result  is  necessary  for  the 
floating/stacking  technique  to  be  broadly  applicable.  Note  as  well  that  the  average  enhancement 
of  peak  area  is  some  13x  for  the  floating/stacking  technique  relative  to  “standard”  injection 
methods. 
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Table  2.9:  Comparison  of  areas  of  electrophoretic  peaks  resulting  from  pinch-&-pullbackand 
floating/stacking  injection  and  separation  of  12  different  eTag™  reporters. 


Peak  AREA 

Peak  AREA 

ratio 

eTag 

Pinch  &  pullback 

Floating/stacking 

P&PB/FS 

1 

15.2 

210 

13.8 

2 

3.1 

43.7 

13.9 

3 

2.8 

37.0 

13.3 

4 

3.3 

41.7 

12.8 

5 

5.2 

72.6 

13.9 

6 

11.4 

145 

12.6 

7 

2.4 

37.3 

15.5 

8 

5.4 

72.7 

13.4 

9 

18.0 

213 

11.8 

10 

8.6 

106 

12.3 

11 

14.9 

169 

11.3 

12 

11.5 

127 

11.0 

average 

1+ 

c/> 

D 

13.0±1.3 

Implementation  of  “Point  Concentrator”  Features  for  Improved  Limits  of  Detection  in 
Electrophoretic  Analysis 


Improvement  of  the  limits  of  detection  in  electrophoretic  separations  offers  the  potential  to  allow 
more  accurate  medical  diagnostics,  forensic  measurements,  and  other  bioassays  with  ever 
smaller  sample  quantities. 

Electrophoresis  is  commonly  used  to  move  ionic  species  for  subsequent  separation  and  detection. 
A  common  problem  with  this  technique  is  that  the  concentration  of  species  at  the  detector  is 
often  too  small  to  detect  with  good  signal-to-noise  ratios.  One  previously  explored  method  to 
address  this  problem  is  the  technique  of  “stacking”  (preconcentrating)  sample  species  by 
changing  the  conductivity  of  the  background  electrolyte  relative  to  that  of  the  sample.  ACLARA 
developed  a  new  technique  introducing  local  geometrical  features  (asperities)  in  the  channels 
used  to  carry  out  electrophoresis,  distorting  the  electric  field  in  the  injector’s  offset  region. 
These  geometrical  changes  caused  a  local  concentration  buildup  of  the  sample  species  to  be 
detected.  The  geometrical  features  explored  include  the  placement  of  triangular  asperities  in  the 
offset  region.  The  overall  geometry  and  placement  of  electrodes  used  in  the  electrophoresis 
model  is  shown  in  Figures  2.51  and  2.52. 
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Figure  2.51:  Channel  and  reservoir  configuration  utilized  for  point  concentrator  studies. 


Figure  2.52:  Channel  geometry/layout  for  point-concentrator  simulations,  showing  the  20  x  31  pm 
triangular  asperity  used  to  create  local  distortion  of  the  electric  field  during  the  loading/injection  process. 

The  sample  channels  were  coated  with  1%  PEO;  the  running  buffer  was  25  mM  HEPES  at  pH 
7.38;  the  sample,  placed  in  the  well  at  electrode  V4,  consisted  of  1  pM  fluorescein  in  25  mM 
HEPES  +  25  mM  NaCl.  A  geometrical  configuration  that  was  found  to  substantially  increase  the 
local  electric  field  strength  in  the  offset  region  is  shown  in  Figure  2.52. 


Simulation  of  the  electrophoresis,  consisting  of  injection  followed  by  separation,  was  made  using 
this  geometry,  and  an  electropherogram  (Figure  2.53)  was  generated.  This  electropherogram  was 
compared  to  an  electropherogram  obtained  by  using  the  same  injection/separation  scheme,  but 
without  a  triangular  step-shaped  asperity  in  the  channel.  The  comparisons  were  done  with  use  of 
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the  floating  scheme  described  in  Table  2.10,  as  well  as  with  a  pinch-and-pullback  configuration 
described  in  Table  2.11. 

Table  2.10:  Timing/voltage  conditions  for  floating/stacking  operation  of  channel  with  triangular 

asperity  in  the  region  of  the  injector. 

V1  y2  yz  V4  Time  (sec) 

Inj  float  500  float  0  12 

Sep  0  float  700  float  25 


Table  2. 1 1 :  Timing/voltage  conditions  for  pinch-and-pullback  operation  of  channel  with 
triangular  asperity  in  the  region  of  the  injector. 


V1 

V2 

V3 

V4 

Time  (sec) 

Inj 

0 

500 

0 

0 

12 

Sep 

0 

380 

700 

380 

25 

The  electropherogram 

shows 

that  the  triangular 

asperity  (“triangle”) 

gives  about 

a  2x 

improvement  in  peak  height  relative  to  the  result  without  the  triangle  in  floating/stacking  mode 
(float),  and  a  much  higher  peak  than  the  use  of  the  triangular  offset  in  pinch-&-pullback  mode. 
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Figure  2.53:  Comparison  of  single-species  electropherogram  using  the  geometrical  configuration  shown 
in  Figure  5.3  (“triangle”)  with  a  floating/stacking  voltage/timing  configuration;  the  same  voltage/timing 
without  the  triangular  asperity  (“fioaf ’);  and  the  use  of  the  asperity  with  the  “pinch-and-pullback” 

configuration. 
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2.4  Biokinetic  Data  Extraction  for  Array-Based  Biosensors  (University  of  Utah  and  CFD 
Research  Corporation) 

Interpreting  binding  kinetics  from  array-based  biosensors  requires  new  data  analysis  tools.  The 
methods  must  be  fast  in  order  to  interpret  the  large  amounts  of  reaction  data  and  accurate  enough 
to  account  for  transport  effects  that  occur  because  of  analyte  concentration  gradients  created 
within  large  flow  cells.  In  order  to  develop  the  appropriate  analysis  tools  for  array-based 
biosensors,  the  team  focused  the  initial  work  on  validating  data  extraction  methods  used  for 
single  flow  cell  systems.  Three-dimensional  CAD  within  CFD-ACE-i-  was  used  to  generate  a 
realistic  construction  of  an  operating  BIACORE  2000  biosensor  flow  cell  (see  Figure  2.54). 
Binding  data  were  then  simulated  under  a  variety  of  conditions,  including  increasing  the 
association  rate  for  the  reaction  from  10,000  to  30,000,000  M'^s'^  (see  Figure  2.54). 


Figure  2.54:  Left  panel  shows  CFD-ACE+  simulation  of  three-dimensional  flow,  diffusion,  and  binding 

reaction  within  a  single  flow  cell. 


Right  panel  shows  response  data  simulated  for  reactions  with  increasing  association  rates  {ka  =  le4,  le5, 
3e6  and  3e7  M"'s'')  but  with  the  same  dissociation  rate  {kd  =  0.03  s'*).  The  red  lines  represent  a  fit  to  the 
data  using  a  fast  two-compartment  model  to  describe  transport  and  binding. 


These  simulated  data  were  then  fit  using  a  fast  two-compartment  model  that  describes  transport 
of  an  analyte  to  the  surface  and  a  reversible  binding  reaction  with  an  immobilized  ligand. 
Excellent  fits  were  obtained  to  the  data  simulated  with  the  different  association  rate  constants  (as 
shown  by  an  overlay  of  the  red  lines  with  the  simulated  data  in  Figure  2.54).  More  importantly, 
the  association  rate  constants  that  were  determined  from  the  data  extraction  program  correlated 
very  well  with  the  values  use  in  the  simulation.  Furthermore,  the  analysis  method  also  returned 
accurate  values  for  the  dissociation  rate  constant  for  all  of  the  simulations  (ka  =  ~0.03  s  ').  The 
analysis  program,  which  utilizes  the  simple  two-compartment  transport  model,  runs  very  quickly 
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and  reduces  analysis  time  by  ~  1000-fold  for  a  comparable  data  set  that  is  analyzed  using  the  full 
three-dimensional  fitting  routine  developed  for  this  work  using  ACE-I-.  Increased  speed  will  be 
particularly  important  when  it  becomes  common  to  analyze  hundreds  of  data  sets  obtained  from 
array-based  biosensors. 

In  order  to  establish  validated  test  reagents  for  molecule  devices,  the  team  characterized  a 
number  of  molecule  interaction  systems  ranging  from  antibody/antigen,  drug/target  and 
DNA/DNA  interactions.  Examples  of  data  from  some  of  the  molecule  systems  the  team 
characterized  are  illustrated  in  Figure  2.55.  This  work  is  critical  for  several  reasons.  Before 
initiating  these  studies,  if  someone  wanted  to  test  the  performance  of  a  device,  there  were  no 
standardized  systems  to  choose  from.  The  team  now  has  sets  of  commercially  available  reagents 
which  can  be  used  by  the  research  community  for  systems  evaluation  and  training.  The  team  has 
published  several  papers  describing  the  characterization  of  some  of  these  systems  and  have 
utilized  them  for  training  and  to  conduct  four  global  benchmark  studies  [58-63].  Having 
standardized  systems  also  allowed  us  to  develop  and  validate  a  novel  approach  to  kinetic  analysis 


Figure  2.55;  Example  data  sets  for  six  validated  molecular  systems. 


(A)  DNSA/CAII,  (B)  BODIPY/CAII,  (C)  GFP/mAb,  (D)  FITC/mAb,  (E)  Protein  A/mAb  and  (F)  a 
DNA/DNA  interaction.  Black  lines  represent  the  responses  collected  from  the  biosensor.  Red  lines 
represent  the  global  fit  of  each  data  set  to  a  simple  1:1  interaction  model. 


that  is  called  “kinetic  titration” [64].  The  data  shown  in  Figure  2.56  are  an  example  of  a  kinetic 
titration  series  carried  out  on  an  antibody/antigen  system  that  had  been  validated  previously.  In  a 


75 


kinetic  titration  assay  analyte  samples  are  tested  from  low  to  high  concentration  without 
regenerating  the  surface.  Kinetic  titrations  allow  one  to  collect  the  required  concentration 
dependent  binding  information  faster  than  traditional  kinetic  methods  which  require 
regenerations  between  binding  cycles.  Kinetic  titrations  will  be  particularly  important  as  array  - 
based  biosensor  technology  is  further  developed,  by  allowing  the  user  to  analyze  a  wide  range  of 
kinetic  systems  without  the  need  for  regeneration.  This  new  method  also  required  the 
development  new  data  fitting  methods  which  have  made  available  to  the  research  community 
through  CFDRC’s  distribution  of  new  software  (Figure  2.56). 

The  research  produced  key  innovations  in  the  areas  of  electrokinetic  stirring  and  electrothermal 
flow  quantitation  and  physics.  The  group  demonstrated  that  one  can  achieve  an  order  of 
magnitude  improvement  in  heterogeneous  reactions  using  AC  electrokinetic  microstirring.  It 
also  demonstrated  that  micro-PlV  can  be  used  to  quantify  ac  electrokinetic  flow,  and  that  use  of 
the  team’s  2  particle  strategy  can  significantly  reduce  errors  due  to  dielectrophoresis. 


Figure  2.56:  Left  panel  shows  kinetic  titration  data  for  PSA/mAb  interaction. 

Right  panel  shows  modeling  program  interface  which  was  developed  and  used  to  fit  titration  data  [64]. 


In  addition  to  the  team’s  own  interest  in  developing  validated  systems,  the  team’s  work  also 
received  the  attention  of  the  National  Institutes  of  Standards  and  Technology.  Members  of  NIST 
saw  a  need  for  developing  standardized  reagents  for  molecular  devices  and  the  team  worked 
closely  with  them  to  design,  develop  and  validate  a  new  system  which  they  could  mass  produce 
and  distribute.  A  bivalent  mutant  of  the  protein  Bamase  was  designed  which  would  produce  a 
system  with  two  binding  sites  containing  different  affinities  (Figure  2.57).  Having  one  system 
with  two  affinities  makes  this  standardized  system  useful  for  testing  a  wider  range  of  devices. 
The  group  is  currently  in  the  final  stages  of  testing  these  reagents  but  have  already  used  them  to 
develop  a  novel  method  of  assaying  long  association  phases  on  Biacore  biosensors  [65]. 
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Barstar 


Figure  2.57;  Concentration  dependent  binding  profiles  for  the  NIST  bivalent  Bamase/Barstar  system. 


Finally,  the  team’s  standardized  reagents  are  being  used  by  a  number  of  instrument  manufactures 
to  benchmark  their  own  technology.  Companies  such  as  Lumera  (Bothell,  WA),  BioRad 
(Hercules,  CA),  and  Biacore  Inc.  (Piscataway  NJ)  requested  and  received  reagents  from  the  team 
for  their  developmental  work.  BioRad  and  BIACORE  later  adopted  two  of  the  systems  produced 
by  the  team  as  training  and  test  reagents  that  they  distribute  with  their  instruments.  The 
experiments  are  currently  using  a  validated  antigen/mAb  system  to  test  the  sensitivity  of  a  new 
array  based  biosensor  from  Lumera.  Therefore,  the  current  efforts  in  developing  and  validating 
molecular  reagents  have  impacted  a  wide  area  of  research  and  will  continue  to  do  so  in  the 
future. 

2.5  Electrokinetic  Instability  (Stanford  University) 

The  team  discovered,  formulated,  and  created  models  of  a  new  electrokinetic  instability.  This 
instability  is  in  fact  an  electrohydrodynamic  instability  that  occurrs  in  the  regime  applicable  to 
electrokinetic  microsystems  and  which  couples  with  electroosmotic  flow  to  produce  convective 
motion.  This  section  desribes  some  background  on  this  work,  presents  sample  experiments, 
introduces  the  current  formulation  of  this  problem,  and  presents  modeling  results  and  their 
comparison  with  experiments. 


2.5.1  Importance  of  Electrokinetic  Instability  Flow  Studies 

Over  the  past  decade  there  has  been  extensive  research  into  the  design  of  devices  that  perform 
chemical  analyses  in  micro-fabricated  fluidic  channel  structures.  Often  referred  to  as  micro  total 
analysis  systems  (//TAS),  these  systems  exhibit  a  mass  transport  regime  that  is  often  distinct 
from  that  of  macro-scale  flow  devices.  Many  of  these  devices  apply  electrokinetic  liquid-phase, 
bioanalytical  techniques  including  capillary  electrophoresis  and  isoelectric  focusing,  and  often 
manipulate  samples  having  poorly  characterized  or  unknown  electrical  conductivities.  As  a 
result,  conductivity  mismatches  often  occur  between  the  sample/reagent  species  and  the 
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background  electrolyte.  In  the  presence  of  applied  electric  fields,  conductivity  gradients  can 
induce  electrohydrodynamic  coupling,  which  can  in  turn  generate  complex,  unstable  flow-fields. 
Flows  exhibiting  these  physics  have  been  reported  in  the  classical  electrohydrodynamics 
literature  (see  for  example,  the  seminal  paper  by  Hoburg  and  Melcher,  1976  [66],  and  a  later 
work  by  Baygents  and  Baldessari,  1998  [67]).  In  this  section,  flow  instabilities  due  to  electric 
field  and  conductivity  gradients  coupling  in  electrokinetic  systems  are  reviewed.  Electrokinetic 
flows  are  a  subset  of  electrohydrodynamics  characterized  by  the  presence  of  an  electrical  double 
layer  and  regimes  where  transport  due  to  molecular  diffusion  is  important. 

Although  desirable  for  rapid-mixing  applications,  electrokinetic  instabilities  are  unwanted  in 
microfluidics  applications  such  as  sample  injection,  separation,  and  controlled  diffusion-limited 
reaction  processes  where  minimal  sample  dispersion  is  required.  This  motivates  research  toward 
a  better  understanding  of  the  conditions  necessary  for  the  onset  of  electrokinetic  flow  instability. 


2.5.2  Brief  Review  of  Work  Conducted  as  Part  of  the  Current  Project 


In  2001  Oddy  et  al.  (Stanford)  [8]  first  reported  observation  of  electrokinetic  instability  (EKI)  in 
a  microchannel  system.  These  experiments  were  performed  in  4-mm-long  glass  capillaries  with 
rectangular  cross  sections,  and  the  instabilities  were  in  general  of  temporal  nature  [8].  In  a 
slightly  different  geometry  (microfluidic  T-junction),  Chen  and  Santiago  also  reported  spatial 
amplification  of  disturbances  which  was  later  identified  as  convective  instability  [68,  69].  In  all 
of  these  experiments,  conductivity  gradients  were  in  the  spanwise  direction  (perpendicular  to  the 
electric  field),  and  there  existed  critical  values  of  the  applied  streamwise  electric  field  above 
which  instabilities  and  significant  stirring  occurred. 

Following  these  initial  experimental  observations,  there  has  been  a  development  of  models  for 
electrokinetic  flow  instabilities.  Models  are  useful  in  predicting  threshold  conditions  for 
instability  onset  as  well  as  other  flow  features  including  coherent  wave  structures  and  mixing 
rate.  Lin  et  al.  [70]  and  Chen  et  al.  [69]  showed  that  a  generalized  EHD  modeling  framework 
(derived  from  the  so-called  “leaky  dielectric”  model  first  developed  by  Melcher  and  Taylor, 
1969  [71])  can  be  used  to  describe  both  low-conductivity,  non-diffuse  charge  dynamics  of 
classical  EHD,  and  the  more  recently  reported  flow  instabilities  of  high-conductivity  electrolyte 
in  electrokinetic  microsystems.  Lin  et  al.  [70]  analyzed  the  temporal  stability  properties  based  on 
a  two-ion  model,  comprising  the  conductivity  transport  equation  along  with  the  conservation  of 
momentum  and  electromigration  current.  They  showed  that  the  model  provided  good  qualitative 
and  fair  quantitative  agreement  with  regard  to  the  threshold  electric  fields  and  critical 
wavenumbers.  Lin  et  al.  [70]  also  presented  non-linear  simulations  of  their  set  of  governing 
equations  that  capture  the  high  Peclet  (or  the  so-called  electric  Rayleigh)  number  stirring  events 
observed  in  experiments.  Using  a  convective  framework,  Chen  et  al.  [69]  showed  that  EKI  could 
manifest  itself  convectively  in  the  presence  of  a  strong  electroosmotic  flow.  In  the  latter  analysis, 
EKI  is  modeled  using  a  linearized,  thin-layer  limit  of  the  Navier-Stokes  equations  coupled  with 
conservation  equations  for  electrical  conductivity  and  current.  The  model  reveals  both 
convectively  and  absolutely  unstable  eigenmodes.  More  recently,  Storey  et  al.  [72]  presented  a 
depth-averaged  version  of  the  governing  equations  used  by  the  Lin  et  al.  model[70].  Their  depth- 
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averaged  model  compared  favorably  with  a  complete  three-dimensional  model  for  thin  channel 
geometries. 


In  the  next  two  sections,  salient  results  from  these  various  studies  are  presented.  A  brief 
presentation  of  sample  experimental  results  is  followed  by  a  depth-averaged  model  useful  in 
predicting  flow  instabilities  and  computing  dispersion  in  electrokinetic  microfluidic  systems. 


2.5.3  Example  Flow  Visualizations  and  Experimental  Results  for  Electrokinetic  Flow 
Instabilities 


Here  are  shown  a  few  exemplary  results  from  these  experiments.  The  microchannel  consisted  of 
a  borosilicate  glass  capillary  (Wilmad-Labglass,  NJ)  with  a  rectangular  cross-section.  The 
channel  was  4  mm  long  (in  the  x  or  streamwise  direction),  1  mm  wide  (in  the  y  or  spanwise 
direction),  and  0. 1  mm  deep  (in  the  z  or  depth  direction).  Two  buffer  streams  initially  occupied 
the  upper  and  lower  halves  of  the  microchannel  resulting  in  a  diffuse  conductivity  gradient  along 
the  spanwise,  y  -direction;  an  electric  field  was  subsequently  applied  along  the  streamwise,  x  - 
direction.  The  conductivity  of  the  buffer  streams  were  50  and  5  //S/cm,  respectively,  resulting  in 
a  conductivity  ratio  of  y  =  l0.  For  flow  visualization,  an  electrically  neutral,  heavy-molecular- 
weight  dye  (VOkDalton)  composed  of  a  dextran-rhodamine  B  conjugate  (Molecular  Probes,  OR) 
was  added  to  the  high  conductivity  buffer  stream.  The  imposed  electric  potential  initiated  an 
electroosmotic  flow  in  the  channel  and,  for  electric  fields  above  a  threshold  value,  electrokinetic 
instabilities. 

Here  are  presented  a  few  exemplary  results  from  the  current  experiments.  A  schematic  of  the 
microchannel  system  used  in  many  of  the  current  parametric  studies  is  shown  in  Figure  2.58 
below.  The  microchannel  consisted  of  a  borosilicate  glass  capillary  (Wilmad-Labglass,  NJ)  with 
a  rectangular  cross-section.  The  channel  was  4  mm  long  (in  the  x  or  streamwise  direction),  1 
mm  wide  (in  the  y  or  spanwise  direction),  and  0.1  mm  deep  (in  the  z  or  depth  direction).  Two 
buffer  streams  initially  occupied  the  upper  and  lower  halves  of  the  microchannel  resulting  in  a 
diffuse  conductivity  gradient  along  the  spanwise,  y  -direction;  an  electric  field  was  subsequently 
applied  along  the  streamwise,  x -direction.  The  conductivity  of  the  buffer  streams  were  50  and  5 
p  S/cm,  respectively,  resulting  in  a  conductivity  ratio  of  y  =  10  .  For  flow  visualization,  an 
electrically  neutral,  heavy-molecular- weight  dye  (VOkDalton)  composed  of  a  dextran-rhodamine 
B  conjugate  (Molecular  Probes,  OR)  was  added  to  the  high  conductivity  buffer  stream.  The 
imposed  electric  potential  initiated  an  electroosmotic  flow  in  the  channel  and,  for  electric  fields 
above  a  threshold  value,  electrokinetic  instabilities. 
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Figure  2.58:  Schematic  of  the  microchannel  setup. 

Two  buffer  solutions  with  differing  electrical  conductivities  are  advected  through  the  microchannel  using 
a  syringe  pump,  resulting  in  a  single  buffer  stream  having  a  spanwise  electrical  conductivity  gradient. 
Images  were  captured  halfway  between  the  inlet  and  outlet  of  the  capillary.  Fluid  motions  were  observed 
using  an  inverted,  epi-fluorescent  microscope  (Nikon  TE300),  4x  (NA  0.2)  objective  and  a  0.6x 
demagnifying  lens,  resulting  in  an  overall  magnification  of  2.4x.  A  CCD  camera  (CoolSnap  fx,  Roper 
Scientific  Inc.,  AZ)  with  a  12-bit  intensity  digitization  resolution  recorded  the  images.  Image  signal-to- 
noise  ratio  SNR  was  improved  by  binning  individual  CCD  pixels  forming  4x4  super  pixels  having  final 

dimensions  of  26.8  x  26.8  pm  in  the  image  plane. 


This  experimental  setup  allowed  for  rapid  ehange  of  parameters  including  the  applied  electric 
field;  the  characteristic  width  of  the  diffuse  region  between  the  two  conductivity  streams  (an 
important  initial  condition);  the  location  of  imaging;  and  the  ratio  of  conductivity  of  the  two 
streams.  The  setup  also  allowed  for  rapid  change  between  rectangular  capillaries  of  a  variety  of 
sizes  (with  cross-section,  inner  dimensions  ranging  from  10  um  to  2  mm)  and  aspect  ratios 
(varying  from  1:1  to  10:1). 

A  representative  set  of  images  from  experiments  conducted  at  1,2,  and  3  kV  applied  potentials 
are  shown  in  Figure  2.59.  The  potentials  were  applied  over  a  distance  of  40  mm,  such  that  they 
were  equivalent  to  applied  fields  of  25,000,  50,000  and  75,000  V/m,  respectively.  In  each  case, 
the  top  figure  of  each  series  shows  the  initial,  undisturbed  interface  between  the  dyed  and  undyed 
buffer  streams  in  the  channel  (t  =  0).  The  successive  images  in  each  column  show  the  temporal 
evolution  of  the  imaged  dye  under  a  constant,  DC  potential.  In  this  color  scheme,  blue 
corresponds  to  the  undyed,  low-conductivity  stream,  and  red  to  the  dyed  high  conductivity 
stream.  For  an  applied  field  of  25,000  V/m,  the  interface  was  only  slightly  perturbed  and  only 
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slight  fluctuations  are  apparent  in  the  images  captured  at  4.0  s  and  5.0  s.  At  the  two  higher 
applied  voltages,  the  interface  exhibited  a  rapidly-growing  wave  pattern  within  the  first  1  s.  The 
unstable  fluid  motion  in  the  channel  buckled  the  interface  and  proceeded  to  stretch  and  fold 
material  lines.  The  transverse  and  fluctuating  velocities  associated  with  this  unstable  motion 
resulted  in  rapid  mixing  of  the  two  streams.  At  the  75,000  V/m  applied  field,  the  channel  reached 
a  well-stirred  state  and  with  nearly-homogeneous  concentration  fields  observable  within  5  s. 

A  representative  set  of  images  from  experiments  conducted  at  1,2,  and  3  kV  applied  potentials 
are  shown  in  Figure  2.59.  The  potentials  were  applied  over  a  distance  of  40  mm,  such  that  they 
were  equivalent  to  applied  fields  of  25,000,  50,000  and  75,000  V/m,  respectively.  In  each  case, 
the  top  figure  of  each  series  shows  the  initial,  undisturbed  interface  between  the  dyed  and  undyed 
buffer  streams  in  the  channel  ( t  =  0 ).  The  successive  images  in  each  column  show  the  temporal 
evolution  of  the  imaged  dye  under  a  constant,  DC  potential.  In  this  color  scheme,  blue 
corresponds  to  the  undyed,  low-conductivity  stream,  and  red  to  the  dyed  high  conductivity 
stream.  For  an  applied  field  of  25,000  V/m,  the  interface  was  only  slightly  perturbed  and  only 
slight  fluctuations  are  apparent  in  the  images  captured  at  4.0  s  and  5.0  s.  At  the  two  higher 
applied  voltages,  the  interface  exhibited  a  rapidly-growing  wave  pattern  within  the  first  1  s.  The 
unstable  fluid  motion  in  the  channel  buckled  the  interface  and  proceeded  to  stretch  and  fold 
material  lines.  The  transverse  and  fluctuating  velocities  associated  with  this  unstable  motion 
resulted  in  rapid  mixing  of  the  two  streams.  At  the  75,000  V/m  applied  field,  the  channel  reached 
a  well-stirred  state  and  with  nearly-homogeneous  concentration  fields  observable  within  5  s. 
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Figure  2.59:  Sample  images  from  the  experiment,  shown  for  applied  fields  of  25,000,  50,000,  and  75,000 
V/m,  corresponding  to  the  first,  second,  and  third  column. 


Images  obtained  at  various  times  are  shown  for  each  column.  The  electric  field  and  bulk  flow  directions 
were  from  left  to  right.  High  voltage  was  applied  as  a  Heaviside  function  at  f  =  0  s.  Each  image 
corresponds  to  a  physical  area  1  mm  wide  (y)  and  3.6  mm  long  ( x ).  The  depth  of  the  channel  is  100  jum 
along  the  z  -direction  (into  the  page).  Small  amplitude  waves  at  /  =  1  s  grew  and  led  to  rapid  stirring  of  the 
initially  distinct  buffer  streams.  The  instability  stretched  and  folded  material  lines  and,  after  about  4  s  for 
the  75,000  V/m  applied  field,  resulted  in  a  well-stirred,  relatively  homogeneous  dye  concentration  field. 

The  time  of  the  images  in  each  row  are  shown  in  the  figure. 

The  team  demonstrated  the  electrokinetic  instability  in  a  variety  of  flow  geomoetries  including 
T-channels  and  X-intersctions  (see  Figure  2.60).  An  example  result  from  across-charmel 
electrokinetic  flow  injection  geometry  is  shown  in  Figure  2.60.  The  figure  shows  experiments  in 
a  50  um  wide  (20  micron  deep)  borosilicate  microchip  system  from  Micralyne  (Alberta, 
Canada).  On  the  left,  is  a  sample  image  showing  a  stable  flow  when  the  two  streams  have  the 
same  conductivity.  On  the  right,  are  two  sample  instantaneous  images  for  the  case  where  the 
center  stream  (from  the  west  to  east  wells)  conductivity  is  100  times  that  of  the  background 
sheath  stream  from  the  north  and  south  wells. 
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Figure  2.60;  Electrokinetic  instability  during  electrokinetic  “pinching”  at  an  X-intersection. 

(a)  Schematic  and  concentration  field  in  the  conductivity-matched  condition,  (b)  Unstable,  fluctuating 
concentration  field  for  a  sample  conductivity  100  times  smaller  than  the  background  electrolyte.  Same 
DC  voltage  scheme  but  conductivity-gradient  case  results  in  rapid  (0.1  sec  scale)  fluctuations  in  a 
region  within  a  few  channel  widths  from  the  intersection.  A  180  V/cm  field  with  ReL  <0.1  is  shown 
(smallest  structures  are  <  5  pm). 


It  is  also  important  to  mention  that,  early  on  in  the  EKI  research  effort,  the  team  designed  and 
fabricated  an  EKI-based  micromixer  which  successfully  demonstrated  rapid  mixing  for  low 
Reynolds  number  flows.  [1]  The  mixer  design  won  the  2001  inventor’s  award  from  the  National 
Inventors  Hall  of  Fame  (established  by  the  U.S.  Patent  Office),  and  a  patent  is  pending.  The 
device  represents  an  enabling  technology  for  binding  and  surface-hybridization  assays  of  “Lab- 
on-a-Chip”  systems. 


2.5.4  Electrokinetic  Flow  Instability  Model  and  Sample  Results 

A  general  theoretical  formulation  developed  in  Lin  et  al.  [70]  is  presented  here.  Using  these 
equations  results  from  various  linear  analyses  as  well  as  nonlinear  simulations,  and  assess  their 
qualitative  and  quantitative  agreements  with  experimental  data  are  shown.  As  a  conclusion  of 
this  section  of  this  technical  report,  the  Stanford  group’s  latest  development  in  a  depth-averaged 
framework  suitable  for  the  study  of  generalized  electrokinetic  flows  in  microchannels  with  thin 
channel  geometry  is  introduced. 

Below  is  a  summary  of  the  governing  equations  and  the  key  parameters  of  interest  in  the 
Stanford  studies.  The  equations  of  the  Stanford  model  are  derived  from  the  conservation  laws  for 
a  dilute,  two-species  electrolyte  solution  [73],  and  (with  justification)  a  relaxation  assumption  to 
simplify  the  equations  is  adopted.  The  scaling  analysis  and  derivations  are  discussed  in  detail  in 
[70]  and  should  not  be  repeated  here.  The  (dimensionless)  equations  read 


^  +  vVcr  =  — VV, 
dt  Ra^ 


(2.49) 
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V-(oVO)  =  0, 


(2.50) 


V^O  =  -p^,  (2.51) 

V  •  V  =  0,  (2.52) 


Re 


\ 


d\ 

—  +  V  Vv 

ydt  j 


=  -Vp  +  V^v  - 


(2.53) 


where  a  is  the  conductivity,  v  is  the  bulk  fluid  velocity,  O  is  the  electric  field  (which  includes 
both  applied  and  generated  components),  and  p  is  pressure.  The  electric  Rayleigh  number 
(similar  to  a  Peclet  number)  and  the  Reynolds  number  that  arise  from  the  nondimensionalization 
are  defined  as: 


(2.54) 


Ra„  = 


D  ’ 


Re  = 


UH 


V 


Here  H  is  half  channel  width,  D  is  an  effective  diffusivity  of  the  conductivity,  and  v  is  the 
kinematic  viscosity  of  the  electrolyte  solution  (usually  aqueous).  The  velocity  the  so-called 
electroviscous  velocity,  is  a  velocity  scale  derived  by  setting  equal  the  electric  body  force  and  the 
viscous  force  in  the  momentum  equation  [66,  69,  70] (Note:  More  discussions  on  these 
parameters  can  also  be  found  in  these  literatures.) . 

First,  a  two-dimensional  formulation  is  presented,  where  it  is  assumed  that  the  flow  exists  only  in 
the  x-y  plane,  and  that  there  are  no  dynamics  in  the  z  -direction.  This  analysis  will  capture  the 
basic  physics  of  the  instability  mechanisms  due  to  the  conductivity  gradient.  A  linear  stability 
analysis  is  used  to  predict  the  regimes  where  it  would  be  expected  rapid  mixing  to  occur.  The 
base  states  are  a  diffused,  spanwise  conductivity  profile  (y)  with  a  (high-to-low) 

conductivity  ratio  of  10,  and  a  shear  electroosmotic  flow  u^=u^(y).  Note  that  the  spanwise 
dependence  of  the  later  was  induced  by  that  of  the  former  via  the  dependence  of  zeta  potential  on 
conductivity.  An  assumption  is  that  the  disturbance  is  periodic  in  the  x  direction,  and  the 
growth  of  amplitude  is  exponential  in  time.  The  results  are,  for  each  streamwise  wave  number  k 
and  applied  field  E^,  a  set  of  eigenvalues  (the  exponential  growth  rates),  together  with  their 
respective  eigenfunctions.  Figure  2.61  is  a  contour  plot  of  the  growth  rates  of  most  unstable 
eigenfunction  in  the  wave  number-Rayleigh  number  (electric  field)  parameter  space.  The  symbol 
5  denotes  the  real  and  dimensional  growth  rate.  The  neutral  stability  curve  is  obtained  by  setting 
j  =  0.  A  threshold  electric  field  is  successfully  captured  from  the  minimal  value  of  on  the 
neutral  stability  curve.  As  originally  reported  by  Baygents  and  Baldessari  [67],  the  team  found 
that  the  inclusion  of  the  diffusive  term  V^alRa^  in  equation  (2.49)  is  crucial  for  the  existence  of 
the  neutral  stability  curve. 
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Figure  2.61;  Contour  plot  of  growth  rates  (s)  versus  wave  number  and  Rayleigh  number. 

Dimensional  applied  electric  field  is  provided  on  the  right  axis.  For  the  case  plotted  here,  the  ratio  of  the 

conductivity  between  the  two  streams  is  10. 


In  addition  to  the  two-dimensional  linear  stability  analysis  presented  above,  the  group  also 
solved  the  full  nonlinear  governing  equations  numerically  to  capture  the  nonlinear  evolution  of 
the  instability  observed  in  the  experiments.  The  initial  conditions  are  the  base  states  plus  a  small 
white  noise  perturbation.  The  solution  details  are  documented  in  Lin  et  al.  [70]. 

Figure  2.62  shows  the  nonlinear  evolution  of  the  simulated  dye  at  various  instances  in  time  and 
for  three  different  electric  fields.  The  model  reproduces  many  of  the  essential  features  observed 
in  the  experiments,  including  the  shape  and  initial  break-up  dynamics  of  the  interface,  the 
transverse  growth  of  a  wave  pattern  in  the  interface,  and  the  roll-up  of  scalar  structures  observed 
at  later  times.  Note  the  similarity  in  the  most  unstable  (and  most  apparent)  wave  number  at  later 
times  between  the  simulation  and  experiments. 


85 


Figure  2.62:  Snapshots  of  the  simulated  dye  field  at  various  instanees  in  time  for  different  driving  electric 

fields. 


In  this  color  scheme  red  corresponds  to  the  high  conductivity  buffer,  and  blue  to  the  low  one.  Each 
column  indicates  a  different  applied  field  and  the  rows  within  each  column  present  selected  snapshots  in 
time.  The  image  corresponds  to  a  physical  domain  of  3.6mmx  1mm.  The  time  for  noticeable  waves  to 

develop  is  decreased  as  the  field  is  increased. 

Despite  similarities  between  wave  number  and  the  dynamics  of  the  interface  breakup,  the 
threshold  imposed  fields  from  both  the  linear  and  nonlinear  predictions  are  lower  than  those 
shown  for  the  experiment  in  Figure  2.59.  For  example,  compare  the  evolution  of  the  dye  at 
25,000  V/m  from  the  experiments  (Figure  2.59,  column  1)  and  the  simulation  (Figure  2.62, 
column  3).  The  results  show  that  the  simulation  at  25,000  V/m  predicts  a  well-stirred  flow  field 
in  less  than  three  seconds  while  experiments  show  that  the  flow  is  stable  on  the  time-scale  of  the 
experiments.  The  simulation  of  25,000  V/m  is  qualitatively  similar  to  the  experimental  flow  at 
75,000  V/m  (Figure  2.59,  column  3).  Despite  the  discrepancy  in  the  magnitude  of  the  applied 
field,  the  Stanford  group  simulation  captures  a  threshold  field  and  scalar  features  qualitatively 
similar  to  the  experiment.  In  the  following  section  possible  causes  for  the  under-prediction  of  the 
threshold  electric  field  is  addressed  by  including  three-dimensional  effects. 


In  comparison  with  the  temporal  instability  analysis  presented  here  (which  is  consistent  with  the 
experimentally  observed  instability  in  the  previous  section),  Chen  et  al.  [69]  analyzed  the 
instability  in  a  convective  framework  which  is  consistent  with  the  spatial  growth  that  was 
observed  in  T-junction  microfluidic  channels.  Among  other  contributions,  the  authors  found  that 
an  important  dimensionless  group  ,  defined  as  the  ratio  of  the  electroviscous  to  electroosmotic 
velocity,  is  critical  in  demarcating  the  absolute/convective  instability  boundary.  Interested 
readers  are  referred  to  [69]. 
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The  previous  sections  have  been  providing  a  two-dimensional  framework  which  appears  to 
capture  the  primary  physics  of  the  current  flow.  However,  the  primary  flaw  in  that  model  is  the 
two-dimensional  assumption,  for  a  channel  with  an  aspect  ratio  of  <5  =  dIH  =  0.1 ,  where  d  denotes 
the  channel  half-width.  Such  a  thin  channel  geometry  was  also  used  in  the  experimental  work  of 
Chen  et  al.[68,  69].  In  three  dimensions,  an  EDL  forms  not  only  on  the  top  and  bottom  walls  of 
y  =  ±\,  but  also  along  the  side  walls  (z  =  +1 ),  and  strongly  drives  the  flow  due  to  the  small  depth 
of  the  channel.  The  three  dimensional  nature  of  a  thin  channel  has  the  added  dynamics  that  the 
fluid  motion  in  the  interior  of  the  channel  is  directly  coupled  to  the  top  and  bottom  wall 
boundary  condition  (determined  in  part  by  the  local  value  of  ion  density). 

In  Lin  et  al.  [70]  the  group  presented  3D  linear  analysis  and  found  that  the  predicted  threshold 
field  was  an  order-of-magnitude  higher  than  that  from  the  2D  linear  analysis,  in  much  closer 
agreement  with  the  experimental  observations.  However,  a  depth-averaging  approach  is 
preferred  since  in  general  the  3D  analysis  (simulations)  are  computationally  more  expensive,  and 
considering  that  the  microchannels  of  Stanford’s  interest  are  “shallow”  in  the  depth  (  z  ) 
direction.  In  Chen  et  al.[69]  a  depth-averaged  analysis  was  performed  on  a  set  of  linearized 
governing  equations,  and  the  resulted  linear  equation  system  was  used  for  convective  instability 
analysis.  In  Storey  et  al.  [72]  a  more  complete  Hele-Shaw  type  of  integrated  momentum  equation 
was  used.  Both  models  yielded  favorable  quantitative  results  when  compared  with  experimental 
data. 

Here  the  ideas  that  were  first  developed  by  Chen  et  al.  [69]  and  Storey  et  al.  [72]  are  extended 
and  completed.  A  generalized,  nonlinear  depth-averaged  model  suitable  for  the  study  of 
electrokinetic  microchannel  flows  with  thin  channel  geometries  is  developed.  It  is  accomplished 
this  through  a  complete  asymptotic  analysis  of  the  full  three-dimensional  equations  based  on  a 
smallness  parameter  which  is  the  channel  cross-sectional  aspect  ratio  S  .  Stanford’s  general 
methodology  follows  a  combined  lubrication  (for  the  momentum  equations)  and  Taylor-Aris  (for 
the  convective-diffusion  of  the  conductivity  field)  approach.  Without  presenting  the  details  of  the 
derivation,  the  final  equations  are: 
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Here  the  overbar  denotes  a  depth-averaged  quantity,  the  operator  denotes  the  in-plane 
gradient  (to  distinguish  from  the  full  three-dimensional  gradient),  and  U  =  u  -  is  the 
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difference  between  the  total  depth-averaged  velocity  and  the  electroosmotic  velocity.  The  main 
contributions  of  this  new  equation  set  are  the  Taylor-dispersion-type  term  in  the  conductivity 
equation,  and  the  Darcy-Brinkman-Forchheimer  (DBF)  type  of  momentum  equation  which  is  of 
second-order  consistency  in  S  [69,  74],  Preliminary  results  in  the  assessment  of  the  validity  and 
accuracy  of  the  model  are  presented. 

Figure  2.63  compares  the  linear  stability  results  for  growth  rate  of  disturbances  versus  wave 
number  at  a  single  Rayleigh  number  of  Ra^  =  5, 000 .  Linear  analyses  is  performed  using  the 
following  three  models:  1.  The  linear  integrated  momentum  equations  used  by  Storey  et  al.  [72]; 
2.  The  linearized  three-dimensional  equations [70];  and  3.  The  depth-averaged  DBF  formulation 
presented  here.  Note  that  at  the  linear  regime  the  Taylor  dispersion  in  Equation  (2.55)  drops  as  a 
higher-order  term,  and  the  only  difference  between  the  models  are  within  the  momentum 
equations.  It  is  found  that  all  three  models  are  in  good  agreement  for  wave  numbers  below  about 
4.  However,  when  compared  with  the  more  accurate  three-dimensional  analysis,  the  DBF 
momentum  equation  provides  significantly  better  results  at  higher  wave  numbers  than  the  lower- 
order,  integrated  momentum  approximation. 


Figure  2.63:  Comparison  of  growth  rates  of  disturbances  as  predicted  by  three  models. 

Shown  here  are  the  real  part  of  the  growth  rates  versus  the  wave  number  for  Ra^  =  5,000  as  computed  with 
the  DBF  momentum  equation  presented  here,  the  integrated  momentum  equation  ( approximation 
[72]),  and  the  three-dimensional  equation  set.  The  DBF  formulation  represents  in-plane  viscous  stresses 
that  quench  unphysical  high  wave  number  growth  and  is  in  agreement  with  the  three-dimensional  result. 


The  characteristics  of  the  DBF  momentum  equations  also  make  it  more  advantageous  to  use  in 
nonlinear  simulations  when  compared  with  the  integrated  momentum  equation.  In  particular,  the 

inclusion  of  in-plane  diffusion  V^u  preserves  a  mathematical  structure  similar  to  the  original 
Navier-Stokes  equations,  and  enables  reproduction  of  boundary  effects  (e.g.,  at  y  =  ±l  walls) 
that  are  not  captured  by  lower-order  approximations.  This  will  be  discussed  further  in  future 
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work,  and  here  are  simply  shown  some  sample  results  of  the  full  nonlinear,  depth-averaged 
simulations  with  Equations  (2.55-2.58).  That  is,  a  model  with  the  combined  effects  of  Taylor 
dispersion  and  the  DBF  momentum  equation.  The  team  tries  to  reproduce  the  experimental 
image  presented  in  Figure  2.59  at  the  two  lower  voltages  (25,000  and  50,000  V/m),  the  result  is 
shown  in  Figure  2.64.  Again  the  model  reproduces  essential  features  observed  in  the  experiments 
such  as  fastest  growing  wave  numbers  and  the  growth  rates  of  the  interface  disturbance 
amplitude.  However,  note  that  the  computations  are  now  at  exactly  the  same  field  strength  as 
those  applied  in  the  experiments  (as  opposed  to  the  unnaturally  low  fields  used  for  comparison 
with  the  simple  2D  model  results  of  Figure  2.62).  Future  work  will  also  include  the  application 
of  the  model  to  different  flow  configurations  such  as  those  used  in  field  amplified  sample 
stacking  (FASS). 


50,000  V/m 


t  =  5  s 


Figure  2.64:  Nonlinear  simulation  of  the  depth-averaged  equation  system. 


This  model  includes  the  combined  effects  of  the  Taylor  dispersion  and  the  DBF  momentum  equation.  In 
the  strongly  nonlinear  regime  the  Taylor  dispersion  acts  as  an  extra  smoothing  mechanism.  The  results 
compare  favorably  with  the  experimental  data  presented  in  Figure  2.62. 


As  shown  in  Figure  2.64  above,  the  depth-averaged,  non-linear  model  of  EKI  quantitatively  the 
threshold  field  for  the  instability  as  a  function  of  conductivity  ratio.  The  model  also  qualitatively 
captures  the  development  of  the  concentration  field  and  dispersion  dynamics. 


2.5.5  Summary  of  Electrokinetic  Instability  Research 

Experimental,  numerical,  and  analytical  results  that  explain  the  basic  mechanisms  behind  an 
electrokinetic  mixing  phenomenon  observed  in  microfluidic  channels  have  been  presented;  this 
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was  a  major  result  and  accomplishment  of  this  project.  The  group  developed  analysis  and 
computations  based  on  various  sets  of  assumptions  for  electrokinetic  flows  in  a  long,  thin 
channel  with  a  transverse  conductivity  gradient.  The  models  developed  here  are  able  to  predict 
general  trends  in  the  data,  as  well  as  many  of  the  qualitative  and  quantitative  aspects  of  the 
observed  flow  field. 

The  models  presented  in  this  work  are  useful  in  optimization  studies,  as  parameter  space  can  be 
spanned  in  simulations  more  quickly  than  in  the  laboratory.  Work  described  by  Oddy  et  al.  [75] 
has  demonstrated  that  oscillatory  electric  field  can  potentially  drive  even  more  vigorous  mixing. 
The  models  presented  here  can  be  used  to  optimize  the  form  of  the  forcing  function,  to  design  the 
shape  of  a  micromixer,  and  to  develop  optimal  control  strategies  for  both  micro-mixing  and  the 
suppression  of  instabilities. 

This  work  has  had  a  significant  impact  on  the  microfluidics  communities.  For  example,  the 
Stanford  group’s  first  three  papers  in  this  area  [69,  70,  76]  have  already  been  referenced  by  other 
journal  articles  over  150  times. 

2.6  Slip  Flow  Investigation  (UCSB) 

In  addition  to  binding  enhancement,  a  significant  amount  of  time  was  spent  investigating  slip 
flow.  The  team  developed  a  simple  two-phase  continuum  model  which  predicts  the  Navier 
boundary  condition  as  a  function  of  thickness  of  gaseous  bubbles  formed  at  a  microchannel 
surface.  The  results  predict  a  slip  length  of  0.9  microns  for  a  18  nano-meter  thick  bubble  layer. 
The  resulting  model  is  also  in  agreement  with  measurements  of  Zhu  &  Granick  [34],  and 
consistent  with  atomic  force  microscope  measurements  of  bubble  formation  by  Tyrell  and 
Attard[l].  The  group  is  devising  experiments  to  test  the  validity  of  the  simple  experimental 
model,  and  to  develop  more  complicated  models  describing  the  non-homogeneous  characteristics 
of  the  velocity  flow  field. 

Prior  measurements  with  micron-resolution  particle  image  velocimetry  (p-PIV)  have  shown  an 
apparent  fluid  slip  for  de-ionized  water  flowing  through  30  x  300  micron  channels  that  are 
coated  with  hydrophobic  octadecyltrichlorosilane  (OTS).  Figure  2.65  compares  the  experimental 
velocity  profiles  for  flow  near  the  wall  of  a  clean  hydrophilic  microchannel  (squares)  and  for 
flow  near  the  wall  of  a  hydrophobic  microchannel  (triangles).  For  a  hydrophobic  channel,  the 
results  show  fluid  slip  near  the  wall  of  approximately  8.5%  of  the  free  stream  velocity. 
Assuming  Navier’ s  hypothesis  (vsiip=P(dv/dy)waii)  as  the  appropriate  boundary  condition  at  the 
wall,  the  measured  slip  corresponds  to  a  slip  length,  P,  of  approximately  0.9pm.  The  mechanism 
for  the  generation  of  apparent  fluid  slip  is  unknown.  Ruckenstein  and  Rajora  (1983)  suggest  that 
fluid  slip  may  develop  from  entrained  and  soluble  gases  forming  a  gap  near  the  wall.  Recently, 
Tyrell  and  Attard  [1]  imaged,  with  an  atomic  force  microscope,  a  hydrophobic  glass  surface 
submerged  in  water.  Their  results  showed  the  presence  of  pancake  shaped,  20  to  30  nm  thick 
nano-bubbles  completely  covering  the  surface.  In  addition,  they  showed  that  the  hydrohobic 
surface  acts  as  a  nucleation  site  pulling  dissolved  gasses  out  of  solution.  Within  10  to  20  minutes 
after  scraping  the  surface  clean,  the  surface  was  once  again  completely  covered  with  nano¬ 
bubbles.  To  determine  if  dissolved  gasses  could  generate  the  apparent  fluid  slip  in  Figure  2.65, 


90 


flow  between  to  infinite  parallel  plates  assuming  nano-bubbles  form  an  effective  air  gap  near  the 
wall  is  examined  analytically.  From  this  simple  one  dimensional  model  one  can  calculate  the 
bubble  height  or  air  gap  required  to  generate  measured  slip  lengths  assuming  Navier’s  hypothesis 
effectively  describes  a  thin  air  gap  at  the  surface. 


For  two  phase  flow  between  two  infinite  parallel  plates  with  a  thin  air  gap  of  thickness,  5,  at  the 
wall  and  a  fluid  layer  of  thickness  2h,  the  solution  for  the  velocity  in  the  water  phase  assuming 
Stokes  flow  for  both  phases,  no  stress  at  the  centerline,  continuity  of  stress  and  velocity  at  the 
air-water  interface,  and  no-slip  at  the  air-wall  interface  is: 
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where  pa  and  pw  are  the  viscosity  of  air  and  water  respectively,  (-dp/dx)  is  the  pressure  drop  and 
y=0  at  the  air  water  interface.  Assuming  Navier’s  hypothesis  effectively  describes  the  thin  air 
gap  at  the  wall  the  velocity  is  set  equal  to  the  slip  length  times  the  shear  rate  and  obtain  an 
equation  for  the  air  gap  thickness  (or  bubble  height)  as  a  function  of  plate  separation  and  slip 
length  : 
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Figure  2.65:  Velocity  profiles  for  water  flow  through  a  hydrophilic  (squares,  no-slip)  and  a  hydrophobic 

(triangles,  slip)  microchannel. 
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Figure  2.66  shows  the  air  gap  thickness  required  for  a  given  plate  sepration  to  generate  a  given 
slip  length.  As  the  plate  separation  decreases  the  slip  length  generated  by  a  given  air  gap 
thickness  increases.  At  very  small  separations  the  slip  length  is  strongly  dependent  on  the  air 
gap  thickness.  In  this  work  is  calculated  a  slip  length  of  0.92  pm.  With  30pm  channels,  this 
yields  a  required  air  gap  thickness  of  approximately  1 8nm.  This  thickness  is  consistent  with  the 
measured  bubble  heights  of  Tyrrell  and  Attard[l].  In  addition,  when  surfaces  are  hydrophobic, 
Zhu  and  Granick  [34]  calculated  slip  lengths  of  approximately  2.5pm  for  a  20nm  separation. 
The  corresponding  bubble  height  of  ~20nm  is  consistent  with  the  required  bubble  height  in  this 
work  and  the  measured  values  of  Tyrrell  and  Attard  [1].  Pit  et.  al.  [77]  calculate  a  slip  length  of 
approximately  400nm  for  a  190pm  separation,  which  would  require  a  bubble  height  of  9nm. 
These  results  suggest  that  the  measured  apparent  fluid  slip  may  develop  as  a  result  of  dissolved 
or  entrained  gasses  accumulating  on  hydrophobic  surfaces.  However,  the  assumption  that  the  air 
gap/nano-bubble  may  be  modeled  as  a  continuum  is  quite  coarse.  At  a  lengthscale  less  than 
approximately  lOOnm,  the  gas  is  rarefied  and  slip  would  occur  at  the  gas/wall  interface.  From 
kinetic  theory  of  gasses  the  velocity  at  the  wall  would  be: 
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where  a  is  the  accommodation  coefficient  and  X  is  the  mean  free  path  of  the  gas.  Solving  the 
two  phase  flow  problem  described  above  but  with  the  kinetic  theory  boundary  condition  yields 
the  following  equation  for  the  slip  length  as  a  function  of  air  gap  height  and  plate  separation: 
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where  s  =Z(4-2a)/(3a).  For  the  experimental  conditions  of  Tretheway  and  Meinhart  [33],  Zhu 
and  Granick  [34],  or  Pit  et.  al.  [77],  Equation  (2.62)  produces  slip  lengths  substantially  larger 
than  those  observed  experimentally.  In  fact,  the  experimentally  measured  slip  lengths  are  not 
obtainable  with  Equation  (2.62).  Thus,  the  assumption  that  the  nano-bubbles  completely  cover 
the  surface  and  form  a  thin  air  gap  along  the  wall  is  incorrect.  However,  if  the  bubbles  partially 
cover  the  surface,  the  wall  would  contain  regions  of  slip  and  no-slip.  Since  the  nano-bubble 
dimensions  of  Tyrell  and  Attard  [1]  are  substantially  smaller  than  the  current  measurement  area, 
the  measured  velocities  average  the  regions  of  slip  and  no-slip.  Thus,  the  appropriate  equation 
for  the  slip  length  as  a  function  of  bubble  height  and  plate  separation  is: 
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here  x  is  the  surface  fraction  covered  with  bubbles.  From  this  equation  is  calculated  a  required 
surface  fraction  covered  with  bubbles  equal  to  .25  to  produced  the  measured  slip  length  of  0.9pm 
for  flow  through  a  microchannel.  This  low  surface  coverage  may  lead  to  substantial  regions  of 
no-slip.  Further  examination  of  the  velocity  profdes  combined  to  form  Figure  2.67  verifies  this 
idea.  Figure  2.67  shows  the  velocity  profiles  along  the  length  of  the  channel.  In  different 
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regions,  the  velocity  profiles  show  varying  degrees  of  slip.  This  may  indicate  regions  of  high 
and  low  bubble  surface  concentration. 


plate 


Figure  2.66:  Solutions  of  Eq.  (2),  showing  air  gap  thickness  as  a  function  of  plate  separation  to  produce  a 

given  slip  length. 


Measured  slip  lengths  for  various  researchers  are  indicated.  Tyrrell  and  Attard[l]  measured  bubble 

heights  of  20-30  nm. 
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Streamwise  Position  (|jm) 


Figure  2.61 :  Velocity  profiles  along  the  microchannel  wall.  The  various  profiles  show  regions  of  slip  and 

no-slip. 

2.7  On-chip  Two-Dimensional  Separations  (Stanford  University) 

The  team  developed  a  novel  on-chip  assay  devices  that  combines  isoelectric  focusing  and  zone 
electrophoresis  to  achieve  a  2D  assay  in  l/30th  of  the  time  of  a  traditional  system  (and  with 
10,000x  less  sample).  Based  on  conventional  two-dimensional  electrophoresis,  this  system 
sequentially  combined  isoelectric  focusing  and  capillary  electrophoresis  in  a  single  planar 
microfluidic  device.  Isoelectric  focusing  (lEF)  separates  chemical  species  based  on  their 
isoelectric  points.  At  completion,  sample  species  reach  a  steady  state  distribution  of 
concentrated  bands  within  a  separation  channel.  Due  to  the  focusing,  steady  state  nature  of  lEF, 
the  assay  is  inherently  insensitive  to  dispersion  caused  by  initial  sample  injection.  Capillary 
electrophoresis  (CE)  separates  species  based  on  their  electrophoretic  mobility,  a  characteristic 
that  is  approximately  proportional  to  the  charge-to-mass  ratio  of  the  species. 

The  planar  geometries  used  for  much  of  this  work  were  designed  in-house  and  fabricated  by 
ACLARA  Biosciences  (Mountain  View,  CA)  in  poly(methyl  methacrylate)  (PMMA)  using  an 
imprinting  and  laminating  technique  similar  to  that  reported  in  the  literature  [78].  The  channel 
cross-sections  were  D-shaped  and  measured  200  pm  wide  and  20  pm  deep. 

The  schematic  of  the  microchip  and  assay  procedure  is  depicted  schematically  in  Figure  2.68. 
The  first  dimension  (horizontal  channel  in  the  figure,  from  A  to  C)  is  2.8  cm  long  and  the  second 
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dimension  (vertical  channel  B-W)  is  2.5  cm  long.  Custom  fixturing  was  designed  and  fabricated 
in-house  to  provide  wells  for  inserting  electrodes  and  allow  for  pressure  fdling  of  liquid  into  the 
channels  [79]. 
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Figure  2.68  Separation  algorithm  for  the  2D  lEF-CE  separation. 

The  device  geometry  is  shown  with  arrows  indicating  the  bulk  fluid  motion  during  each  serial  separation 
step.  The  first  dimension  (lEF)  extends  from  reservoir  A  (anolyte)  to  reservoir  C  (catholyte).  The  second 
dimension  (CE)  extends  from  reservoir  B  (buffer)  to  reservoir  W  (waste).  The  imaged  area  is  indicated 
by  the  dashed  box,  D.  The  voltage  applied  during  each  step  of  the  2D  separation  is  shown  in  the  table 

(HV:  high  voltage,  G:  ground,  F:  float). 

Figure  2.68  and  the  CCD  images  presented  in  Figure  2.69  aid  in  description  of  the  2D  separation 
algorithm.  After  the  initial  focusing  step,  the  2D  separation  algorithm  consisted  of  multiple 
iterations  through  the  following  sequence  of  steps:  (Figure  2.69  A)  an  IFF  separation  with 
simultaneous  EOF  mobilization  of  sample  species  along  the  axis  of  the  first  dimension 
(lEF/EOF),  (Figure  2.69  B)  electrokinetic  sampling  of  a  relatively  small  volume  of  the  lEF 
dimension  constituents  into  the  CE  channel,  and  (Figure  2.69  C)  a  CE  separation  in  the  second 
dimension.  Figure  2.69  D  shows  the  start  of  the  next  separation  cycle  in  which  lEF  and 
simultaneous  mobilization  have  been  reinitiated  in  the  first  dimension.  As  this  is  a  serially 
implemented  2D  separation,  the  separation  sequence  is  repeated  until  all  fluid  volumes  from  the 
lEF  dimension  have  been  sampled  into  the  CE  dimension. 
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Figure  2.69:  CCD  images  during  speeies  sampling. 

(A)  Species  are  focused  by  lEF  in  the  first  dimension  (dark  bands  in  horizontal  channel).  Simultaneously, 
the  bands  are  mobilized  towards  the  catholyte  reservoir  by  low-dispersion  EOF.  (B)  Once  a  fluid  volume 
of  interest,  n,  reaches  the  microchannel  intersection,  all  electrodes  are  switched  to  electrically  float  for  3  s. 
(C)  High  voltage  is  then  applied  at  reservoir  B  and  reservoir  W  is  grounded,  initiating  sample  separation 
in  the  second  dimension.  (D)  Upon  completion  of  the  CE  separation,  lEF/EOF  is  reinitiated  causing 
sample  species  to  refocus  and  the  next  fluidic  volume  (n-1)  to  migrate  to  the  intersection.  This  sequence 
is  repeated  until  all  fluidic  volumes  are  sampled  from  first  dimension  into  the  second. 

Simultaneous  focusing  and  EOF-driven  mobilization  transported  focused  sample  species  to  the 
sampling  intersection  with  an  electroosmotic  velocity  of  order  20  pm/sec  (Figure  2.69  A,  D). 
After  a  fluidic  volume  of  interest  reached  the  intersection  (Figure  2.69  B),  all  reservoirs  were 
switched  to  electrically  float  for  3  s  prior  to  sampling.  During  this  period,  focused  sample 
species  began  to  defocus,  as  presumably  did  the  much  lower  molecular  weight  ampholytes.  The 
floating  step  was  incorporated  into  the  voltage  algorithm  to  electrically  decouple  the  two  assay 
dimensions  and  to  allow  for  a  more  homogenized  pH  field  near  the  intersection.  Mass  continuity 
ensures  that  the  bulk  of  fluid  sampled  from  the  IFF  dimension  into  the  CE  dimension  was 
replaced  with  unfocused  ampholyte  solution  (Figure  2.69  C).  This  algorithm  resulted  in  an  lEF 
separation  that  was  essentially  ‘parked’  during  each  CE  analysis,  as  species  were  refocused  {E  = 
350  V/cm  for  6  -  10  sec)  prior  to  additional  CE  analyses  (Figure  2.69  D).  Each  sampling  event, 
CE  separation,  and  lEF  re-focus  cycle  had  a  period  of  45  s. 

Two-dimensional  “gel-like”  plots  were  constructed  from  a  time-sequence  of  CCD  images 
collected  during  the  2D  separations.  The  gel-like  plots  were  formatted  to  display  inverted  gray 
scale  intensity  information  so  as  to  mimic  a  slab-gel  result  with  dark  regions  corresponding  to 
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high  fluorescence  intensity  zones.  Figure  2.70:  shows  the  results  of  a  2D  separation  using  a  gel¬ 
like  format  for  three  time  steps  during  the  2D  analysis  (At  of  3  s,  5  s,  and  7s).  Each  vertical 
column  (‘lane’)  in  the  gel-like  plot  is  generated  from  a  spatial  electropherogram  that  indicates  the 
intensity  of  fluorescence  along  the  CE  separation  channel.  In  each  gel-like  plot,  the 
electropherograms  correspond  to  equal-duration  CE  separations  (i.e.,  the  CE  data  for  each  2D 
plot  were  collected  at  equal  time  periods  after  the  respective  electrokinetic  injection). 
Accordingly,  the  ordinate  axis  of  the  gel-like  plot  corresponds  to  the  axial  coordinate  along  the 
second  dimension  of  the  assay.  Since  the  method  developed  in  this  work  is  a  serial  2D  analysis, 
the  vertical  lanes  correspond  to  CE  analyses  of  fluid  volumes  sampled  from  different,  adjacent 
locations  along  the  lEF  dimension;  thus,  each  vertical  lane  is  an  analysis  of  a  single  sampling 
event  from  a  discrete  p/ range.  As  a  consequence,  the  abscissa  of  each  2D  plot  is  proportional  to 
the  axial  dimension  along  the  lEF  channel.  Recall  that  during  lEF,  the  focused  bands  are 
mobilized  from  left  to  right  into  the  injection  region  with  a  mobilization  velocity  of  roughly  20 
pm/s.  The  left-most  lane  of  each  gel-like  plot  presented  in  Figure  2.70:  corresponds  to  the  first 
sampling  event,  while  columns  to  the  right  correspond  to  subsequent  sampling  from  fluid 
volumes  containing  decreasing  p/  values. 


Ax  =  3mm  *■ 

lEF  dimension  (mm) 


Figure  2.70:  Gel-like  plots  of  an  lEF-CE  separation  at  CE  analysis  times  of  3  s,  5  s  and  7  s. 


The  horizontal  axis  corresponds  to  the  relative  position  of  each  fluid  element  during  the  lEF  separation. 
Approximately  3  mm,  or  15%,  of  the  total  lEF  channel  length  was  sampled  {E  =  350  V/cm).  The  vertical 
axis  corresponds  to  the  spatial  axial  dimension  of  the  subsequent  CE  separations  {E  =  390  V/cm). 
Species  not  identified  in  companion  ID  separations  are  labeled  as  peaks  1,2,  and  3. 
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2.7.1  2D  Assay  Performance 


The  team  demonstrated  a  single,  planar,  polymer  microdevice  that  serially  integrates  two  rapid, 
orthogonal  chip-based  separation  schemes  (lEF  and  CE).  The  boundary  conditions  of  the 
terminal  channel  reservoirs  (e.g.,  applied  potential,  chemistry,  and  pH)  in  a  cross  geometry  were 
used  to  govern  the  separation  mechanisms  throughout  the  analysis.  A  mixture  of  fluorescent 
sample  species  was  focused  during  lEF  and  mobilized  by  low-dispersion  EOF  to  a  junction, 
where  effluent  volumes  were  repeatedly  electrokinetically  sampled  into  the  CE  dimension.  Such 
electrokinetic  control  of  sample  species  and  bulk  fluid  motion  is  advantageous  for  automated  on- 
chip  systems.  lEF  in  2D  system  showed  rapid  peak  generation  (less  than  1  min),  reduced  EOF 
without  the  introduction  of  high  viscosity  additives,  resultant  highly  concentrated  sample  species 
(70-fold),  and  a  high  peak  capacity.  Species  behavior  in  the  second  dimension  of  the  lEF-CE 
system  was  consistent  with  a  CE  separation  mechanism  (i.e.,  a  constant  velocity  difference 
between  neighboring  peaks  and  diffusional  broadening  of  injected  sample  species).  These  results 
suggest  that  the  two  separation  mechanisms  (lEF  and  CE)  remain  independent  despite  the  fact 
that  the  two  separation  dimensions  are  fluidically  coupled  and  use  the  same  liquid  medium. 
Results  from  the  2D  analysis  suggest  increased  separation  resolution  over  that  of  the 
corresponding  uncoupled  ID  separations,  as  peaks  that  were  not  detectable  during  ID 
separations  under  identical  conditions  were  apparent  in  the  2D  separation.  The  complete  2D 
system  was  estimated  to  have  a  peak  capacity  of  ~1300. 


The  peak  capacity  and  resolution  may  be  further  improved  by  reducing  the  size  of  each  fluid 
volume  sampled  into  the  second  dimension  (perhaps  by  reducing  the  width  of  the  extraction 
junction),  increasing  the  charmel  length  of  the  second  dimension,  and  reducing  dispersion  during 
sample  handling  between  the  first  and  second  dimensions.  The  low  viscosity  of  the  separation 
media,  in  combination  with  the  short  separation  charmel  length  scales  used,  resulted  in  rapid  lEF 
and  CE  analyses.  That  noted,  the  analysis  time  of  this  system  was  dominated  by  that  of  the 
second  dimension  and  could  be  reduced  further  through  implementation  of  a  manifold  of  CE 
charmels  that  would  allow  sampling  and  analysis  of  numerous  analyte  volumes  from  the  first 
dimension  in  parallel.  The  total  analysis  time  for  a  parallel  system  should  decrease  as  the 
number  of  additional  CE  charmels  increases,  although  system  complexity  will  also  increase.  The 
liquid-phase  approach  presented  may  facilitate  further  system  integration,  including  fraction 
collection  and  coupling  of  lEF-CE  with  a  third  dimension  (e.g.,  mass  spectrometry). 


The  combination  of  these  two  powerful  separation  techniques,  in  a  single  system,  can  be  used  to 
obtain  two  important  physical  characteristics  (isoelectric  point  and  mobility),  a  significant 
preconcentration  of  the  sample,  and  a  high  separation  resolution.  These  species  are  well  resolved 
after  less  than  a  minute  of  elapsed  focusing  time.  The  full  characterization  and  performance  of 
the  system  is  described  in  [80]. 
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2. 7.2  Benchmark  Comparisons  with  Competing  Technologies 


Comparisons  between  several  key  performance  parameters  of  the  system  described  here  and  that 
of  various  other  systems  are  presented  in  Table  2.12.  These  benchmark  comparisons  were  made 
to:  classical  2D  slab-gel  electrophoresis,  a  capillary-based  2D  system  involving  lEF,  and  an  on- 
chip  2D  separation  system  involving  CE  as  the  second  dimension. 

Table  2.12:  Multi-dimensional  system  performance  comparison. 

Comparisons  are  made  between  slab-gel  2D  electrophoresis,  a  capillary-based  system  that  couples  lEF 
with  mass  spectrometry,  an  on-chip  system  that  couples  MEKC  with  CE,  and  the  system  developed  in 

this  work. 


2D 

Electrophoresis 

(Slab-gel) 

CIEF-MS 

On-chip 

MEKC-CE 

On-chip 

lEF-CE 

Analysis 

time 

Several  days 

150  min 

10  min 

~60  min 

Volume  of 
first 

dimension 

analyzed 

100% 

— 

10% 

100  % 

1000-3000 

Peak 

capacity 

(standard) 

10-10^ 

>700 

180-360 

-800 

(large  format) 

RSD% 

>  20% 

5-10%  (?) 

— 

>  30%  (?) 

pi  resolution 

0.01 

<0.1  pH  unit 

N/A 

<0.1  pH 
unit 

Required 

sample 

amount 

100  mg  -  >  1  mg 

300  ng 

— 

10  ng 

Automation 

No  completely 
integrated 
systems,  only 
semiautomated 

Possible 

Possible 

Possible 
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2. 7.3  Impact  of  2D  Assay  Work 


The  development  of  robust,  fully-flinctional  multi-dimensional  assay  chips  applicable  to 
complex  analytes  is  considered  by  many  to  be  the  "holy  grail"  in  microfluidic  systems 
development.  The  Stanford  two-dimensional  assay  system  was  the  first  to  combine  capillary 
isoelectric  focusing  (CIEF)  (a  method  of  separating  and  focusing  proteins  in  a  pH  gradient)  and 
CE  on  a  chip  [81],  The  device  was  less  than  2.5  cm  on  each  side  (1.1  mm  thick)  and  automates 
sample  injection  and  separation.  The  initial  lEF  step  focuses  samples  to  more  than  lOOX  their 
concentration,  and  a  pH  range  of  3  to  9  is  addressable.  The  peak  capacity  here  achieved  is 
comparable  to  that  of  commercially-available  macroscale  devices  while  the  assay  time  is  less 
than  1  hr  (vs.  a  typical  2  days  for  the  macroscale  system).  The  required  sample  amount  is  0.05 
micrograms  vs.  order  1  mg  for  the  macroscale  systems.  This  device  and  its  associated 
innovations  have  the  potential  of  replacing  2-D  slab  gel  assays  in  a  variety  of  applications.  The 
system  holds  promise  as  a  basis  for  high-throughput,  high-resolution  protein  and  peptide 
analysis.  The  team  demonstrated  the  performance  of  this  device  in  a  challenging  six  protein 
assay  in  a  collaboration  with  Prof  P.J.  Utz  of  the  Stanford  Medical  School  [81].  The  Stanford 
team’s  close  collaboration  with  Utz  is  impacting  many  areas  of  the  Stanford  team’s  work  and  is 
an  excellent  example  of  the  importance  of  interdisciplinary  research  in  the  field  of  microfluidics. 

2.8  Field  Amplified  Sample  Stacking  Experiments  (Stanford  University) 

The  team  developed  a  novel  PASS  CE  system  with  a  porous  polymer  structure  that  facilitates  the 
establishment  of  high  conductivity  gradients  required  to  achieve  ultra-high  concentration 
increases  using  PASS [4].  The  porous  structure  minimized  instabilities  associated  with  the 
mismatch  in  conductivity.  A  signal  increase  greater  than  a  factor  of  1000  is  obtained,  which,  at 
that  time,  was  the  highest  sensitivity-enhancement  reported  to  date  using  on-chip  PASS  as  a 
stand-alone  method. 

A  porous  polymer  plug  was  fabricated  in  a  glass  microchannel  with  a  double-T  injection 
geometry  (Micralyne,  Alberta,  Canada)  using  photoinitiated  polymerization.  The  porous 
structure  was  photo-defined  at  a  desired  location  using  a  reagent  injection  process  (in 
conjunction  with  a  Mylar  film  mask)  which  the  team  developed. 

The  group  applied  this  PASS  system  to  effect  a  separation  of  2  pM  and  1  pM  initial 
concentrations  of  fluorescein  and  bodipy,  respectively.  Figure  2.71  shows  the  separation  of 
sample  analytes  detected  at  10  mm  downstream  from  the  injection  region.  The 
electropherograms  shown  are  determined  by  spatially  integrating  full-field  detection  data  (over  a 
5000  pm^  region  to  simulate  the  detection  of  a  photo-multiplier  tube  device).  First,  it  is  shown 
that  the  use  of  the  device  in  a  CE  injection  and  separation  with  matched  sample  and  background 
buffers,  i.e.  no  PASS  (Figure  2.72a).  Second,  the  group  have  demonstrated  that  a  signal  increase 
by  more  than  a  factor  of  1000  (Figure  2.72b)  using  an  aggressive  PASS  process  chemistry  with  a 
high  conductivity  (77.6  mS/cm)  background  buffer  and  a  low  conductivity  (60.1  pS/cm) 
background  buffer  and  sample.  Additional  details  on  the  experiments  and  PASS  dynamics  may 
be  found  in  [82]  and  [83]. 
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Boundary  generation  North 


Boundary  generation  North 


South  (a) 


Sample  Loading  North 


Figure  2.71;  Schematic  of  FASS/CE  protocol. 

(a)  High  conductivity  buffer  is  injected  from  the  east  reservoir,  (b)  Low  conductivity 
buffer  is  introduced  from  the  west  reservoir.  Porous  structure  provides  high  fluidic- 
resistance  which  minimizes  the  mixing  of  two  buffers  at  the  boundary,  (c)  Sample  is 
electrokinetically  loaded  in  the  double-T  injector.  Negatively  charged  sample  ions 
electromigrate  from  the  south  reservoir  (grounded)  to  the  north  reservoir  (high  voltage),  (d) 
Sample  stacking,  separation,  and  detection. 
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(a)  (b) 


Figure  2.72:  Electropherograms  of  fluorescein  and  bodipy  separation. 

Fluorescence  signal  was  normalized  with  exposure  time.  The  signal  increase  is  on  the 
order  of  1000  fold  for  the  stacked  case.  The  position  of  detector  is  10  mm  from  the 
downstream  intersection  of  the  channel,  (a)  Electropherogram  of  analytes  in  a  CE 
separation  and  detection  performed  without  stacking.  Conductivity  ratio,  y  =  1,  exposure 
time  is  50  ms.  (b)  Stacked  CE  Electropherogram.  Conductivity  ratio,  y  =  1290,  exposure 
time  is  5  ms.  Difference  in  electromigration  time  is  due  to  the  effects  of  PASS. 
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3.0  MICROFLUIDIC  SYSTEMS  FOR  HIGH  SAMPLE  STACKING  ASSAYS 


3.1  Introduction 

An  additional  task  was  added  to  the  intial  research  effort  that  focused  on  developing  rapid 
sample  pre-concentration  methods  to  improve  the  sensitivity  of  on-chip  assays  (with  a  target 
figure  of  merit  of  10,000  fold  increase  by  project  end).  Leveraging  the  results  on  field  amplified 
sample  stacking  and  electrokinetic  instabilities  from  the  initial  project,  the  Stanford  team  was 
able  to  develop  improved  stacking  methods  using  field  amplified  sample  stacking  (PASS), 
thermal  gradient  focusing,  and  isotachophoresis  (ITP).  Ultimately,  the  team  experimentally 
demonstrated  million-fold  sample  concentration  increase  (approximately  three  orders  of 
magnitude  better  than  the  second  best  ever  demonstrated  on  or  off  chip).  This  section  reviews 
the  results  of  this  experimental  work. 

3.2  Field  Amplified  Sample  Stacking  Results 

As  part  of  the  extension  project,  narrow-channel  chips  were  designed  to  have  a  comparable 
hydraulic  resistance  to  the  porous  plug  chip  design  developed  under  the  original  project[84]. 
Figure  3.1  shows  the  layout  of  the  new  stacking  chip  and  images  of  narrow-channel  regions. 


Figure  3.73:  (a)  Narrow-channel  stacking  chip  geometry. 

Channel  length  between  Rl-Jl,  R2-J1,  and  R4-J1  are  6.5  mm,  length  between  R3-J1  is  45  mm. 
(b)  &  (c)  Reflective  mode  microscopy  images  of  narrow-charmel  with  a  double-T  injection 
geometry.  The  microchannels  are  10  pm  deep  with  a  D-shape  characteristic  of  wet  etching  of 
glass.  Wide-charmels  are  120  pm  wide,  and  narrow-channel  is  40  pm  wide  at  the  top.  The 
lengths  of  narrow-channels  are  200  pm  (b)  or  500  pm  (c). 


High  conductivity  ratio  PASS  experiments  were  performed  using  these  narrow-channel  plug 
microchips.  A  12,000-fold  concentration  increase  was  achieved  through  the  use  of  a  modified 
sample  loading  scheme  that  was  developed  (Figure  3.2).  The  new  sample  loading  scheme 
enabled  the  introduction  of  sufficient  sample  ions  for  ultra-high  conductivity  ratio  sample 
stacking  (Figure  3.3).  Various  EOF  suppression  reagents  and  pretreatment  methods  were  also 
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explored,  and  it  was  determined  that  a  poly-A^-hydroxyethylacrylamide  (PHEA)  coating  achieved 
the  greatest  EOF  reduction.  The  electroosmotic  mobility  of  a  glass  microchannel  was  reduced  to 
2%  of  the  original  value  by  adding  0.1%  PHEA  in  buffer  solution  (10  mM  4-(2-Hydroxy ethyl) 
piperazine- 1-ethanesulfonic  acid,  HEPES,  pH  =  7). 


Figure  3.74:  Plot  of  signal  peaks  for  sample  stacking  of  AlexaFluor  analyte. 


The  highest  peak  shows  a  stacking  ratio  (i.e.,  concentration  increase)  of  12,000-fold.  The  initial  sample 
concentration  was  0.1  pM  of  AlexaFluor  in  DI  water  (a  =  4.72  pS/cm),  and  the  high  conductivity  solution 
was  0.5  M  KCl  in  5  mM  HEPES  (59.1  mS/cm).  The  conductivity  ratio  was  12.5x10^. 


104 


STEP  1. 


low  CJ 
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STEP  3. 


GND 


Figure  3.75:  Schematic  of  PASS  assay  protocol. 


(i)  The  north,  east,  and  west  reservoirs  are  filled  with  low  conductivity  buffer,  high 
conductivity  buffer,  and  sample  solution,  respectively.  The  conductivity  boundary  is  formed 
by  applying  vacuum  to  the  south  reservoir,  (ii)  PASS  is  initiated  by  applying  high  voltage 
and  ground  to  the  east  and  west  reservoirs,  respectively.  The  yellow  arrow  shows  the 
direction  of  electric  field.  Negatively  charged  samples  electromigrate  towards  the  anode 
under  suppressed  EOF.  (iii)  The  voltage  is  switched  to  stop  further  sample  introduction. 


The  12,000-fold  stacking  represents  a  likely  practical  limit  for  stacking  with  PASS.  At  higher 
conductivity  ratios,  both  nonlinear  effects  and  electrokinetic  instabilities  were  observed.  These 
non-linear  effects  showed  behavior  similar  to  that  of  ITP,  another  electrokinetic  pre¬ 
concentration  technique.  The  results  from  the  PASS  studies  suggested  that  ITP  itself  might  yield 
improved  stacking  performance.  The  Stanford  team  pursued  this  approach  and  the  results  are 
summarized  later  in  this  section. 

3.3  Isotachophoretic  Stacking  Results 

Isotachphoresis  is  another  electrokinetic  preconcentration  technique  where  sample  ions  stack  in 
adjent  regions  according  to  their  electrophoretic  mobility,  ultimately  creating  a  train  of  sample 
zones  which  migrate  at  a  constant  speed  (hence  “isotacho”).  Based  upon  the  insights  from  the 
PASS  work,  Stanford  initiated  development  of  an  ultra  high  stacking  ITP/CE  method  (Figure 
3.4). 
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Figure  3.76:  Schematic  of  ITP/CE  assay  protocol. 

Configurations  of  co-ions  are  also  shown  at  each  step,  (i)  The  north  and  the  south  reservoirs  are 
filled  with  leading  electrolyte  (LE),  and  the  west  reservoir  is  filled  with  a  mixture  of  trailing 
electrolyte  (TE)  and  sample.  The  TE/LE  boundary  is  formed  by  applying  vacuum  at  the  south 
reservoir.  White  arrows  show  the  direction  of  pressure-driven  flows,  (ii)  ITP  stacking  is  initiated  by 
applying  high  voltage  and  ground  at  the  east  and  west  reservoirs,  respectively.  The  black  arrow 
denotes  the  direction  of  electric  field.  Sample  ions  (negatively  charged)  electromigrate  toward  the 
anode  as  EOF  is  suppressed.  The  early  stage  of  ITP  stacking  results  in  a  partial  separation  (i.e., 
moving  boundary  electrophoresis),  (iii)  The  field  is  switched  toward  the  north  reservoir  to  initiate 
CE.  ITP  stacking  occurs  until  the  leading  ions  overtake  the  trailing  ions,  (iv)  Separation  of  sample 
analytes  occurs  further  downstream  where  sample  ions  electromigrate  in  homogeneous  background 
electrol3de. _ 


The  hybrid  ITP/CE  method  did  not  require  a  specialized  microchannel  design  or  complex  flow 
control.  It  could  therefore  be  implemented  with  easily  available  “off-the-shelf’  chip  systems 
using  common  voltage  control  methods  and  buffer  chemistries.  High  sample  stacking  was 
achieved  through  a  variety  of  optimizations.  These  included:  using  a  high  concentration  leading 
electrolyte  and  low  initial  sample  concentration,  as  suggested  by  non-dispersive  ID  theory; 
suppressing  EOF  to  minimize  dispersion  due  to  EOF  mismatch;  leveraging  the  high  conductivity 
ratio  and  using  proper  channel  lengths  to  maximize  electric  field  in  the  sample  zone  and  increase 
electric  Peclet  number;  and  implementing  single  interface  stacking  to  yield  a  large  effective 
injected  sample  volume. 

ITP  stacking  is  implemented  using  a  leading  electrolyte  (EE)  with  relatively  high  mobility  ions 
and  a  trailing  electrolyte  (TE)  with  low  mobility  ions.  The  ion  mobilities  of  the  LE  and  TE  are 
respectively  lower  and  higher  than  those  of  the  sample  ions,  so  sample  ions  focus  within  a 
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narrow  zone  between  the  two  electrolytes  and  then  migrate  at  the  same  velocity.  After 
preconcentration,  the  ITP  is  terminated  and  CE  separation  is  initiated  using  a  novel  injection 
method  which  depletes  LE  ions  and  replaces  them  with  TE  ions. 

The  resulting  isotachophoresis  injection  method  was  very  robust  and  easily  integrated  with 
microchip-based  capillary  electrophoresis  (CE)  devices.  The  Stanford  team  was  able  to  improve 
the  performance  of  the  ITP/CE  protocol  and  to  achieve  better  than  one  million-fold  sample 
preconcentration  of  Alexa  Fluor  488  and  performed  high  sensitivity  electrophoretic  separations. 

The  effectiveness  of  the  ITP/CE  protocol  was  demonstrated  with  separations  of  Alexa  Fluor  488 
and  bodipy.  Figure  3.5  shows  the  separations  of  sample  analytes  detected  30  mm  downstream  of 
the  intersection  without  (inset)  and  with  ITP  stacking.  The  electropherograms  were  determined 
by  spatially  integrating  full-field  CCD  imaging  data  over  a  60  by  60  pm  region  centered  on  the 
channel  centerline  to  simulate  the  detection  of  a  point-wise  photodetector.  The  exposure  time 
and  CCD  frame  rate  were  10  ms  and  50  frames  per  second,  respectively.  The  inset  of  Figure  3.5 
shows  the  separation  without  ITP  with  initial  (relatively  high)  concentrations  of  100  nM  and  a 
uniform  background  electrolyte  consisting  of  5  mM  HEPES  buffer  (pH  7.0).  Signal  intensity 
was  normalized  using  the  flatfield  image  signal  of  100  nM  Alexa  Fluor  488.  The  applied  field 
was  280  V/cm.  The  SNRs  of  Alexa  Fluor  488  (first  peak)  and  bodipy  (second  peak)  were  12.9 
and  10.8,  respectively,  and  the  peak  resolution  was  29.4.  SNR  is  defined  as  the  ratio  of  peak 
intensity  to  twice  the  standard  deviation  of  background  noise;  and  resolution  as  the  distance 
between  two  peaks  divided  by  the  full-width  half-maximum  of  the  wider  peak.  Separation  was 
performed  on  a  sample  mixture  diluted  by  a  factor  of  1E5  (1  pM  solutions  each  of  Alexa  Fluor 
488  and  bodipy)  using  the  ITP/CE  method.  The  result  is  shown  in  Figure  3.5.  The  injection/ITP 
time  was  40  s  and  with  a  nominal  applied  field  of  220  V/cm.  The  resolution  of  these  peaks  was 
6. 1  and  the  SNRs  were  respectively  57.6  and  39.0  for  Alexa  Fluor  488  and  bodipy.  The  ITP 
phase  of  the  experiment  shown  in  Figure  3.5  achieved  a  concentration  increase  of  approximately 
0.5  million- fold  immediately  prior  to  initiation  of  the  CZE  mode.  This  experiment  achieved 
concentration  increase  of  6.4E4  relative  to  the  initial  sample  concentration  (1  pM),  and  signal 
increase  of  4.5E5-fold  relative  to  the  unstacked  case. 
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Figure  3.77:  Isotachophoretic  stacking  and  detection  of  Alexa  Fluor  488. 

Figure  3.6  shows  a  spatiotemporal  intensity  plot  of  the  calibrated  fluorescence  intensity  (a 
measure  of  concentration)  versus  separation  distance  and  time  from  which  sample  peak  reaches 
the  detection  window,  for  an  initial  sample  concentration  of  100  fM  of  Alexa  Fluor  488  dye. 

The  LE  was  1  M  NaCl,  and  the  TE  was  5  mM  HEPES  buffer.  The  microscope  is  equipped  with 
a  lOX  objective  (Numerical  Aperture  (N.A.)  of  0.4)  with  viewing  dimensions  of  1.2  by  0.1  mm 
in  the  object  plane.  A  brief  (~  2  min)  sample  stacking  step  enables  the  detection  of  100  fM 
analyte  concentration  with  SNR  =11.  Measured  concentration  distributions  yield  a  maximum 
measured  concentration  increase,  Cl,  of  approximately  2E6-fold.  This  concentration  increase  is 
achieved  within  120  s  and  30  mm  downstream  of  the  injection  point.  To  the  knowledge  of  the 
Stanford  group,  this  is  the  highest  sample  preconcentration  for  either  capillary  or  on-chip 
electrophoresis  systems. 


108 


4 


Figure  3.78:  Detection  of  Alexa  Fluor  488  at  an  initial  concentration  of  100  fM,  greater  than  million-fold 

concentration  increase. 

LE  and  TE  were  respectively  1  M  NaCl  and  5  mM  HEPES.  A  lOX  objective  (N.A.  of  0.4)  was  used. 


3.4  Comparison  of  the  Current  High  Sensitivity  Methods  with  Other  Work 


Here  is  presented  a  comparison  of  the  work  with  the  detection  limits  of  various  point-wise  CE 
detection  modalities  and  associated  research  efforts.  This  comparison  is  an  effort  to  characterize 
and  summarize  the  many  point-wise  detection  methods  and  serves  as  a  short  introduction  to  the 
electroosmotic  flow  diagnostics  systems  discussed  in  this  chapter. 


Figure  3.7  shows  the  detection  limits  of  CE  detection  systems  as  a  function  of  their  molar 
sensitivity.  The  figure  is  adapted  from  a  less  comprehensive  version  published  in  the  book  by 
Landers[85].  Molar  sensitivity  is  plotted  on  the  abscissa,  and  is  formulated  here  as  the  number 
of  sample  moles  loaded  onto  the  capillary  column.  The  ordinate  of  the  graph  shows  the 
concentration  of  the  sample  of  interest.  The  concentration  limit  of  detection  for  a  given  sample 
amount  is  then  dictated  by  the  detection  technique  used.  Trends  are  apparent  for  three  detection 
schemes:  electrochemical,  fluorescence  and  ultraviolet  (UV)  absorbance. 
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Figure  3.79:  Detection  techniques  for  capillary  electrophoresis  (adapted  from  the  CRC  Handbook  of 

Capillary  Electrophoresis  [85],  with  permission). 


Molar  sensitivity  (here  defined  as  the  number  of  moles  which  are  detected)  and  concentration  are  plotted 
on  the  abscissa  and  ordinate,  respectively.  Fluorescence  detection  is  observed  to  have  the  highest 
sensitivity  among  single-point  CE  detection  schemes.  Also  shown  in  the  plot  is  the  detection  capability 
for  fluorescence  detection  in  a  microfiuidic  system  with  an  intensified  12-bit  CCD  camera. 


Electrochemical  detection  includes  amperometry,  conductivity,  and  potentiometry.  The 
drawback  with  this  detection  scheme  is  the  requirement  for  integrated  electrode  structures  within 
the  microchannel.  In  contrast  to  electrochemical  detection,  the  sensitivity  of  UV  absorbance  is 
weak  but  is  still  very  popular  due  to  its  simplicity  of  operation  and  the  fact  that  untagged  samples 
are  detectable.  A  majority  of  the  samples  of  interest,  including  peptides  and  proteins,  absorb  in 
UV  wavelengths.  Despite  its  popularity  in  free-standing  capillary  systems,  UV  detection  is  not 
easily  applicable  to  microchannel  systems  etched  in  planar  glass  substrates  due  to  the  short 
optical  path  lengths  available  from  the  shallow  channel  depths,  and  the  difficulties  associated 
with  transmitted  light  optics  in  an  etched  microchannel.  Example  detection  limits  from  various 
high-sensitivity,  on-chip  detection  efforts  are  shown  as  filled  (blue)  circles. 


no 


Shown  together  with  this  review  of  other  work  are  results  from  the  current  high  sensitivity  ITP 
work.  As  mentioned  above,  the  team  developed  a  high  sensitivity  on-chip  CE  detection  system 
that  leverages  both  optimized  optics  and  with  a  single-interface,  transient  ITP  stacking  protocol. 
The  sensitivity  of  this  system  (without  sample  stacking)  is  shown  by  the  data  point  labeled  “High 
sensitivity  system.”  As  described  by  Jung  et  al.  2006  [84]  high  efficiency  photodetectors  are 
combined,  and  performed  extensive  calibration  studies  of  signal  intensity  dependence  on 
excitation  lasing  power  and  objective  specifications  (N.A.s  and  magnifications).  A  high  N.A. 
objective  was  used  with  relatively  low  magnification  (yet  compatible  with  cover  glass  thickness 
of  microchip)  and  optimization  of  laser  power  so  as  not  to  overly  saturate  sample  analytes. 
Fluorophores  with  high  quantum  yield  and  photo-  and  chemical  stability  were  also  used.  As  an 
example  of  the  sensitivity  of  this  system  the  group  demonstrated  effective  single-molecule 
detection  of  10  pM  Alexa  Fluor  488  under  pressure-flow  conditions  (bulk  velocity  of  1  cm  /s  and 
no  stacking)  using  an  optimized  laser-induced  confocal  fluorescence  microscope  setup  fitted 
with  60x  objective  (N.A.  of  1.4). 

The  sensitivity  achieved  by  combining  more  standard  CCD  type  detection  (i.e.,  the  detection 
system  used  for  Figure  3.5)  with  the  current  optimized  million- fold-capacity  ITP  stacking  is 
indicated  by  the  data  point  labeled  “ITP  stacking  -i-  low  sensitivity  system.”  Lastly,  the 
sensitivity  achieved  by  combining  the  current  high  sensitivity  detection  (i.e.,  the  detection 
system  used  for  Figure  3.6)  with  the  current  optimized  ITP  stacking  is  indicated  by  the  data  point 
labeled  “ITP  stacking  -i-  high  sensitivity  system.”  This  high  sensitivity  CE  detection  system  and 
a  single-interface,  transient  ITP  stacking  method  developed  previously  were  combined  to  show 
separation  and  detection  of  100  aM  concentrations  each  of  Alexa  Fluor  488  and  bodipy.  This  is, 
to  the  team’s  knowledge,  the  highest  sensitivity  capillary  electrophoresis  separation  and 
detection  [84]. 

3.5  Some  Reasons  Why  the  Current  ITP  Stacking  is  Three  or  More  Orders  of  Magnitude 
Better  than  other  ITP 

The  group  developed  a  ID  non-dispersive  model  and  scaling  of  dispersion  analysis[64,  86]. 
Although  these  scaling  arguments  are  simple,  they  yield  important  insight  into  key  ITP  stacking 
parameters  and  suggest  strategies  for  optimizing  ITP  in  practice.  These  strategies  include  using 
high  EE  concentration  and  low  initial  sample  concentration  to  maximize  achievable 
concentration  increase;  suppression  of  EOF  to  minimize  dispersions;  and  implementation  of  a 
single-column  ITP  configuration  (where  initially  there  is  a  single  interface  between  the  EE  and 
the  TE/sample  mixture)  to  inject  a  large  effective  sample  width.  These  strategies  were 
implemented  in  the  design  of  an  ITP/CE  method  that  combines  a  simple  and  robust  single¬ 
column  ITP  stacking  step  with  a  subsequent  CE  step.  First,  a  sample/TE  mixture  was  injected  to 
initiate  ITP  stacking  and  then  injected  LE  at  the  channel  intersection.  The  latter  LE  ion  stream 
overspeeds  the  TE  ions  behind  the  stacked  sample  zone,  terminates  ITP,  and  initiates  CE 
separation. 

The  team  performed  an  experimental  parametric  study  focused  on  the  variation  of  LE,  TE,  and 
initial  sample  concentration.  Consistent  with  the  ID  non-dispersive  model,  it  was  found  that 
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stacked  sample  concentration  is  proportional  to  the  LE  concentration,  and  the  concentration 
increase,  Cl,  (conveniently)  increases  as  the  initial  sample  concentration  decreases.  It  was  also 
found  that  the  stacked  sample  concentration  is  a  strong  function  of  the  TE  concentration  and  the 
initial  sample  concentration.  For  the  current  configuration,  low  conductivity  TE  zones  result  in 
higher  TE  zone  electric  fields  and  therefore  faster  rates  of  stacking.  High  TE  zone  fields  are 
probably  also  correlated  with  high  stacked  zone  fields.  The  hypothesis  is  that  this  finding, 
combined  with  the  large  effective  injection  length  of  the  single-interface  ITP  configuration, 
results  in  large  sample  electric/electrophoretic  Peclet  numbers  and  efficient  stacking  that  is  less 
susceptible  to  dispersion.  Results  suggest  that  comprehensive  multi-species  models  with  coupled 
fluid  flow  (including  perhaps  effects  of  non-uniform  and  dynamic  zeta  potentials),  current 
conservation,  and  convective-diffusion-electromigration  conservation  are  needed  to  fully 
describe  ITP.  Including  of  reaction  buffer  kinetics  may  also  be  important.  The  latter  is  the 
major  motivation  for  the  Stanford  team’s  collaboration  with  the  Sandia  group  whose  work  is 
described  below. 

3.6  Conclusions  for  High  Sample  Stacking  Assays 

Electrokinetic  stacking  techniques  are  a  powerful  tool  for  preconcentration  of  analytes  of 
interest.  The  Stanford  team  successfully  met  and  exceeded  the  target  goal  of  10,000-fold  sample 
stacking,  and  further  developed  techniques  capable  of  delivering  better  than  million- fold  stacking 
in  less  than  two  minutes.  The  team  also  demonstrated  the  highest  sensitivity  ever  demonstrated 
for  an  electrophoresis  separation:  detection  of  100  aM  samples. 


4.0  MODELING  OF  HIGH  SAMPLE  STACKING  ASSAYS 

4.1  Introduction 

This  section  of  the  report  details  work  on  the  simulation  of  sample  stacking  using 
isotachophoresis.  The  simulations  were  performed  with  a  microchannel  flow  model  that  has  a 
detailed  representation  of  the  various  coupled  physical  phenomena  that  are  relevant  to  sample 
stacking  in  flows  with  high  electrical  field  strength  gradients.  The  model  relies  on  an  equilibrium 
buffer  electrolyte  formulation  to  represent  the  changes  in  the  buffer  composition  as  a  function  of 
the  solution  pR.  The  effect  of  free  charges  in  the  solution,  which  are  non-negligible  in  cases  with 
high  gradients  in  the  electrical  field  strength,  is  taken  into  account  in  the  momentum  and  charge 
conservation  equations. 

Extensive  parametric  studies  were  performed  using  isotachophoresis  sample  stacking  methods 
developed  by  Prof  Santiago’s  research  group  at  Stanford  University.  One-dimensional 
simulations  showed  good  qualitative  agreement  with  the  experimental  observations  in  terms  of 
the  early-time  evolution  of  the  dye  in  the  stacking  zone.  The  simulations  also  revealed  that  the 
buffer  pR  can  change  quite  dramatically  in  this  area.  Despite  the  promising  qualitative 
agreement,  there  remains  a  significant  discrepancy  between  the  simulated  and  experimentally 
observed  stacking  ratios.  Qualitative  differences  were  also  observed  in  the  late-time  steady-state 
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behavior.  It  is  surmised  that  more  detailed  two-dimensional  simulations  will  be  required  to 
investigate  these  discrepancies. 


This  task  on  sample  stacking  was  led  by  Stanford  University  with  Juan  Santiago  as  the  principal 
investigator.  Sandia  National  laboratories  (SNL)  and  Johns  Hopkins  University  (JHU)  operated 
under  subcontract  to  Stanford  for  this  project.  The  goals  of  the  SNL  and  JHU  contribution  were 
to  provide  modeling  support  for  the  experimental  work  in  Juan  Santiago’s  group.  Given  the 
particular  focus  of  the  experimental  work  on  ITP  and  the  evident  high  performance  of  this 
sample  stacking  technique,  the  modeling  studies  focused  specifically  on  sample  stacking  in  the 
context  of  ITP. 


4.2  Formulation 

The  simulation  of  isotachophoresis  for  sample  stacking  in  microfluidic  channels  requires  the 
resolution  of  many  coupled  physical  processes.  For  this  purpose,  the  group  leveraged  a 
multiphysics  microchannel  simulation  code  that  was  developed  as  part  of  a  previous,  DARPA- 
funded  project  [87].  This  code  already  contained  the  capabilities  to  simulate  many  of  the 
physical  phenomena  occurring  in  microchannel  flow,  and  was  enhanced  as  part  of  this  work  in 
order  to  handle  challenges  specific  to  isotachophoresis  with  high  gradients  in  the  field  strength. 

This  section  gives  an  overview  of  the  formulation  of  the  resulting  microchannel  flow  model  and 
highlights  the  aspects  of  the  model  that  were  added  as  part  of  this  project. 

4.2.1  Flow  Field 


The  flow  field  is  described  by  the  two-dimensional  unsteady,  constant  density,  Navier-Stokes 
equations  along  with  the  continuity  equation. 


V-u  =  0 

u  •  Vu  =  -Vp-i-  V  •  |^v|^(Vu)-i-(Vu/ J 


In  these  equations,  u  is  the  velocity  vector,  p  is  the  pressure  (normalized  by  density),  v  is  the 
kinematic  viscosity,  pe  is  the  free  charge  density  and  (p  is  the  electric  potential.  The  term  -yO^V  (p 
represents  the  body  force  exerted  by  the  electric  field  on  free  charges  in  the  solution.  This 
physical  phenomenon  was  added  to  the  model  as  part  of  this  work.  For  the  simulation  of 
electroosmotic  flow,  a  slip  wall  velocity  is  used  as  a  boundary  condition,  with  a  magnitude 
predicted  by  the  Helmholtz-Schmoluchowski  equation: 
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(4.2) 


u  =—V^(p 


where  is  the  wall  velocity,  s  is  the  electrical  permittivity,  //  is  the  dynamic  viscosity,  and  C  is 
the  zeta-potential,  which  is  modeled  as  a  function  of  the  local  pR  and  buffer  molarity. 

4.2.2  Species  Transport 


Species  transport  is  governed  by  the  following  equations: 


^  +  v.[c,(»  +  »')]  =  v.(AVc,) 

K  =  -PihFV(p 


(3.3) 


where  Ci  represents  the  species  concentration  and  Z),  is  the  diffusion  coefficient.  The 
electrophoretic  velocity  u®  depends  on  the  species  mobility  its  charge  number  z,  and  the 
Faraday  constant  F.  The  values  of  the  mobility  and  diffusion  coefficient  are  connected  through 
the  Nemst-Einstein  relationship,/).  =  RT /4^  ,  where  R  is  the  characteristic  gas  constant. 


The  species  transport  equation  above  is  used  for  all  species  that  do  not  participate  in  dissociation 
reactions,  such  as  the  dye  (sample)  ions,  D",  and  the  ions  Na^  and  Cf  of  the  fully  dissociated  salt 
NaCl.  However,  the  components  of  weak  acids  or  salts  that  make  up  the  electrolyte  buffer  in 
microchannels  participate  in  very  fast  electrochemical  dissociation  reactions.  Because  these 
dissociation  reactions  are  much  faster  than  other  transport  processes,  they  are  modeled  with  an 
equilibrium  formulation  [88].  Specifically,  for  the  HEPES  buffer  used  in  the  experiments  for  this 
project,  the  formulation  is  as  follows. 

Consider  the  HEPES  buffer  consisting  of  three  components,  H2L^,  HE*,  and  L',  which  participate 
in  the  following  dissociation  reactions: 

H,L^<  >H^+HL- 

HL-<  >H^+L~ 

Further,  define  the  total  buffer  molarity  ^  =  [H2L^  ]  + [HL“  ]  + [L“  ]  .  Since  the  dissociation 
reactions  do  not  change  the  value  of  this  total  buffer  molarity  (as  those  reactions  are  between  the 
three  ions  that  make  up  the  molarity),  this  total  molarity  can  be  advanced  in  time  using  the 
species  transport  equation  that  is  used  for  non-dissociating  species.  Given  the  total  buffer 
molarity  6,  the  relative  contributions  of  each  component  can  then  be  obtained  from: 
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[«!*]  = 

II 

K,K^+K\H^]  +  [H^f 
K,K^  +K\H^]  +  {H^f 


(4.4) 


These  relationships,  however,  first  require  the  eoneentration  of  H^,  to  be  known.  This 
concentration  is  obtained  from  the  charge  conservation  equation,  =  F^z.  c. ,  which  for  the 
solutions  considered  in  this  work  translates  to: 


[Nan-[Cr]-[D-]  +  [H^]- 


+  \oc 


F 


(4.5) 


where  Kw  is  the  water  dissociation  constant.  The  charge  density  pe  is  calculated  from  the 
Maxwell  equation  for  the  electric  potential  V  •(£V^)=  -p^ .  This  charge  density  is  usually  very 

small  and  can  be  neglected.  However,  for  the  simulation  of  sample  stacking  with  large  gradients 
in  the  electrical  field  strength  £"  =  -V<p ,  the  charge  density  term  can  become  of  the  same  order  of 
magnitude  as  some  of  the  other  terms  in  the  charge  conservation  equation. 

4.2.3  Electric  Potential 


The  electric  potential  (p  is  obtained  from  the  current  continuity  equation: 

Spe 


V  •  (cr V (p')=  -F^z,.V  •  {p^  Vc. )+  u  •  + 

c7  =  F^J^zf/3iC. 


dt 


(4.6) 


where  a  represents  the  electrical  conductivity  of  the  solution,  and  the  summations  go  over  all 
species  in  the  solution.  Note  that  for  the  cases  studied  in  this  work,  the  last  two  terms  in  the 

dp 

current  continuity  equation,  u  •  H - ^ ,  are  generally  very  small  and  can  be  left  out. 

dt 


4.2.4  Implementation 


This  section  discusses  the  implementation  of  the  physics  model  that  was  outlined  in  the  previous 
section.  The  isotachophoresis  (ITP)  configuration  in  the  microchannel  provides  for  sample 
stacking  at  the  interface  between  a  leading  electrolyte  (LE)  region  and  a  trailing  electrolyte  (TE) 
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region.  The  present  implementation  ineludes  NaCl  in  the  LE,  and  both  HEPES  and  a  fluoreseent 
dye  in  the  LE,  where  the  dye  is  the  sample  being  stacked. 

Algorithms  were  developed  for  solving  the  governing  equations  in  a  rectangular,  two- 
dimensional,  microchannel  geometry,  with  inflow-outflow  boundaries.  The  construction  built  on 
the  Sandia  group’s  earlier  work  in  [87],  incorporating  necessary  extensions  for  the  additional 
physical  details  in  the  present  context,  particularly:  (1)  the  discontinuous  HEPES  buffer-NaCl 
configuration,  and  (2)  non-electroneutral  effects.  Highlights  of  the  algorithmic  construction  are 
given  in  the  following.  The  presentation  involves  discussions  of  both  two-dimensional  (2D)  and 
one-dimensional  (ID)  constructions.  The  ID  implementation  was  found  to  be  useful  in 
algorithmic  development  and  testing.  It  was  also  used  for  most  of  the  results  reported  further 
below,  given  its  computational  efficiency  and  robustness. 

The  incompressible  Navier-Stokes  equations  are  solved  using  a  projection  scheme,  where  the 
convective  terms  are  integrated  using  order  AdamS-BASHFORTH  (AB3)  time  integration, 
while  the  dissipation  terms  are  solved  using  fractional  stepping  with  2“'^  order  RUNGE-KUTTA 
(RK2)  time  integration.  The  pressure  field  is  solved  for  using  a  fast-fourier  transform  (EFT) 
direct  solver.  The  velocity  field  and  the  electrostatic  body  force  term,  in  addition  to  the  wall  zeta- 
potential  driven  velocity  boundary  condition,  provide  the  two-way  coupling  between  the 
momentum  solution  with  Navier-Stokes  and  the  species  equations. 

The  species  equations  are  solved  using  a  strategy  that  involves  specialized  handling  for  species 
dependent  on  fast  (equilibrated)  reactions,  as  described  above.  The  solution  strategy  is  intimately 
coupled  with  both  the  charge  conservation  constraint  and  the  electrostatic  field  solution  [73]. 
Time  integration  of  these  equations  is  shown  schematically  in  Figure  4.1,  using  an  RK2 
construction.  The  team  has  also  implemented  an  AB3  construction  for  these  same  equations,  not 
shown  in  this  schematic.  It  has  also  implemented  both  centered  difference  and  Godunov  upwind 
discretizations  of  the  convective  terms.  Diffusive  terms  were  discretized  using  centered 
difference  constructions. 


The  top  row  in  Figure  4.1  shows  the  state  variables  associated  with  each  mesh  cell.  These 
variables  include: 

•  Solute  concentrations  Strongly  dissociated  ion  concentrations  are  stored  explicitly  at 
each  grid  cell.  This  includes  the  Alexa-Fluor  488  dye  [D  ],  Sodium  and  Chloride  ions 
[Na  ^]  and  [Cl'].  Hydrogen  [H*]  and  hydroxide  [OH']  concentrations  are  also  stored 
explicitly.  In  contrast,  weakly  dissociated  ions  are  represented  by  a  single  concentration 
of  the  total  molarity  of  all  dissociation  states  in  a  given  cell  {[HL]  in  the  current  case) 
plus  the  a-coefficients  described  above. 

•  Electric  potential  The  electric  potential  <!>  is  computed  at  each  grid  cell.  For  the  ID  case, 
a  direct  numerical  solution  is  implemented.  For  the  2D  case,  a  multigrid  solver  (mg2d) 
was  used. 
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Bulk  Flow  The  bulk  flow  velocity  is  found  from  the  above  Navier-Stokes  solver  for  the 
2D  case.  In  the  ID  case,  it  is  simply  specified,  as  it  is  dictated  by  boundary  conditions  on 
the  domain. 
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Figure  4.80:  The  state  and  equations  used  to  advance  the  simulation  in  time  for  an  RK2  integrator. 


As  Figure  4.1  shows,  the  electric  field  is  derived  from  the  electric  potential  by  a  simple  gradient 
calculation.  The  electric  field  and  bulk  flow  velocity  are  then  used  to  compute  the  rate  of  change 
of  each  transported  solute  concentration.  Three  effects  are  modeled:  electromigration  of  charged 
species,  diffusion,  and  bulk  flow.  Because  water  is  the  solvent  in  which  the  solutes  are 
transported  -  which  effectively  removes  any  upper  limits  on  [H  and  [OH ']  -  and  because  the 
reaction  kinetics  of  hydrogen  and  hydroxide  ions  are  much  faster  than  the  integration  time  scale, 
these  species  are  handled  separately.  Also,  note  that  while  [19]  is  transported,  the  concentration 
coefficients  are  dependent  on  the  pH  of  the  solution  and  must  be  solved  for  in  tandem  with  [H 
and  [OH  ]  by  applying  a  charge  conservation  to  each  cell  of  the  grid.  This  leaves  us  with  a 
quartic  equation  in  terms  of  [H*].  Once  this  equation  is  solved,  the  species  concentrations  can  be 
moved  forward  in  time  with  an  EULER  forward  step.  Derived  quantities  including  the 
conductivity  of  the  medium  and  charge  diffusivity  are  recomputed  and  used  to  update  the  electric 
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potential.  This  is  the  final  calculation  in  the  first  half  of  the  RUNGE-KUTTA  integration  (“Runge 
Kutta  1”  in  Figure  4.1). 

A  corrector  step,  repeating  all  the  steps  in  the  previous  paragraph  using  the  state  of  the  EULER 
forward  step,  is  shown  in  the  “Runge  Kutta  2”  section  in  Figure  4. 1 . 


Figure  4.2  shows  the  experimental  device  schematic  with  the  computational  domain  highlighted. 
For  computational  efficiency,  the  portion  of  the  chamber  where  sample  stacking  occurs  is 
modeled.  The  electric  potential  is  constrained  to  match  the  applied  voltage  difference  with 
Dirichlet  boundary  conditions.  The  concentration  of  each  transported  species  at  the  left  and  right 
boundaries  is  modeled  with  a  fixed-concentration  condition  when  that  species  is  flowing  into  the 
domain  and  as  a  zero-gradient  when  that  species  is  leaving  the  domain.  The  top  and  bottom  grid 
boundaries  have  zero-gradient  boundary  conditions  applied  to  all  species  concentrations  and  to 
the  electric  potential.  Initial  conditions  are  defined  based  on  experimental  conditions  as  well  as 
charge  concentration  and  buffer  equilibrium  constraints. 


Figure  4.81:  The  experimental  apparatus,  with  the  portion  modeled  by  the  simulation  shown  in  green. 


4.2.5  Results 

In  the  following,  a  presentation  of  the  results  of  the  buffer  model  is  followed  by  ID  and  2D  ITP 
results  in  the  microchannel. 

Buffer  Titration  Curve 

As  mentioned  above,  a  HEPES  buffer  was  used  in  the  sample  stacking  experiments.  To  validate 
the  buffer  model  in  the  current  simulation  code,  a  comparison  was  performed  between  the 
experimental  and  simulated  titration  curve  for  a  20  mM  HEPES  buffer.  Figure  4.82  shows  that 
the  simulated  curve  agrees  very  well  with  the  experimental  data. 
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Sample  Stacking  in  ID  at  4  kV/cm 

This  section  gives  an  overview  of  the  results  of  a  ID  simulation  of  sample  stacking  via 
isotachophoresis  (ITP)  for  a  case  with  nominal  field  strength  of  4  kV/cm  in  the  trailing 
electrolyte  (TE).  In  this  case,  the  TE  was  initialized  with  5  mM  HEPES  buffer  and  1  nM  dye  (D' 
).  The  leading  electrolyte  (LE)  was  initialized  with  1  M  Cf.  The  initial  pH  was  set  to  7  and  the 
initial  Na^  concentration  was  calculated  in  order  to  achieve  electroneutrality.  This  case  was 
simulated  over  a  2  cm  long  domain,  with  the  initial  LE  and  TE  zones  each  taking  up  1  cm.  An 
average  field  strength  of  2  kV/cm  was  prescribed  over  this  domain.  However,  as  the  conductivity 
in  the  LE  zone  is  significantly  higher  than  in  the  TE  zone,  almost  the  entire  potential  drop  occurs 
in  the  TE  zone,  resulting  in  a  nominal  field  strength  in  the  TE  zone  of  4  kV/cm.  These  case 
settings  correspond  closely  to  the  experiments  in  which  stacking  ratios  of  up  to  10^  were 
observed. 


This  case  was  simulated  using  4096  grid  points,  with  a  time  step  of  25  ps,  over  a  total  of  200  s  of 
simulated  time.  The  convection  terms  were  calculated  with  the  centered  difference  scheme  and 
time  integration  was  performed  with  AB3.  While  the  current  implementation  allows  for  non¬ 
electroneutral  conditions  in  the  bulk,  as  indicated  above,  parametric  studies  have  been  conducted 
with  and  without  free  charges.  They  showed  a  relatively  minor  role  of  free  charges  in  the  ID 
context.  In  the  results  presented  below,  no  free  charges  were  considered  in  the  bulk 
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(electroneutrality  enforced  everywhere).  Moreover,  with  no  loss  of  generality,  the  bulk  velocity 
was  set  to  zero  over  the  whole  domain.  Figures  4.4  and  4.5  show  the  profiles  of  the  Na"^  and  Cf 
ion  concentrations  as  a  function  of  the  streamwise  coordinate,  x,  at  5  different  points  in  time. 


Both  the  Na^  and  Cl'  ions  clearly  move  downstream  in  a  front.  Unlike  the  Cl'  ion  concentration, 
the  Na^  ion  concentration  develops  a  second  step  in  its  profile  between  the  initial  location  of  the 
interface  and  the  leading  edge  of  the  front.  The  level  of  this  intermediate  step  corresponds  to  the 
amount  of  the  HEPES  buffer  accumulation  in  this  area.  As  shown  in  Figure  4.6,  even  though  the 
total  molarity  of  the  HEPES  buffer  is  only  5  mM  in  the  TE,  the  buffer  species  tend  to  accumulate 
behind  the  leading  ITP  front.  This  is  because  the  buffer  in  this  area  consists  primarily  of 
negatively  charged  ions  (L  ),  as  shown  in  Figure  4.7.  The  mobility  of  the  L'  ions  is  smaller  than 
the  mobility  of  the  CF  ions  in  the  LE  as  well  as  the  mobility  of  the  dye  (D  )  ions.  The  L'  ions 
therefore  accumulate  in  this  region,  and  their  negative  charge  contribution  is  offset  by  the  local 
Na^  ions. 
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Figures  4.8  and  4.9  show  the  evolution  of  the  concentration  and  the  corresponding  pH.  Even 
though  the  solution  started  out  with  a  uniform  pH  of  7,  the  pH  rapidly  increases  to  a  value  of 
over  11.  This  not  only  affects  the  composition  of  the  HEPES  buffer,  it  can  also  influence  the 
fluorescence  of  the  dyes  used  in  the  experiments. 
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The  evolution  of  the  electrieal  eonductivity  field  is  shown  in  Figure  4.10.  As  the  very  mobile  Na^ 
and  Cr  ions  dominate  the  solution,  it  is  not  surprising  that  the  conduetivity  closely  follows  their 
profile,  with  a  double  step  function.  The  electric  field  strength,  however,  does  not  show  this  step 
function,  at  least  not  on  the  scale  of  Figure  4.1 1.  This  is  because  the  electrical  conductivity  is  so 
much  higher  in  the  LE-stacking  zone,  as  compared  to  the  TE  zone,  that  the  potential  drop  (and 
corresponding  field  strength)  in  the  LE-stacking  zone  are  almost  negligible.  The  transport  of 
positive  ions  (by  diffusion  and  electromigration)  into  the  TE  zone  does  result  in  a  slight  increase 
over  time  in  the  conductivity  in  this  area.  This  results  in  a  decrease  of  the  overall  resistance  of 
the  channel,  with  an  (amplitude)  increase  of  the  field  strength  (and  current)  as  a  consequence. 
The  behavior  of  the  negatively  charged  dye  (D  )  concentration  is  shown  in  Figure  4.1 1. 


Figure  4.10;  Evolution  of  the  conductivity. 


Figure  4. 1 1 :  Evolution  of  the  electric  field 
strength. 
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With  an  initial  concentration  of  just  1  nM,  the  dye  concentration  at  t  =  0  s  is  not  visible  on  this 
graph.  However,  as  soon  as  the  electric  field  is  applied,  the  dye  starts  to  stack  behind  the  Cf 
front,  and  after  50  s,  its  concentration  has  already  increased  2000-fold.  After  200  s,  the  dye 
stacking  ratio  is  a  little  larger  than  11,000.  Note  that  even  after  200  s,  the  stacking  continues 
without  any  sign  of  leveling  off. 


The  simulated  ion  concentration  profdes  for  this  case  compare  qualitatively  very  well  with  the 
experimentally  observed  profdes.  Nevertheless,  it  is  clear  that  quantitatively,  the  observed 
stacking  ratio  of  11,000  for  this  case  is  much  lower  than  the  experimentally  observed  ratio  of 
close  to  lO'’.  One  difference  with  the  experimental  conditions  is  that,  in  the  simulation,  the  initial 
concentration  of  the  dye  was  set  to  1  nM,  while  the  experiments  used  an  initial  concentration  of 
1  pM  in  order  to  achieve  such  high  stacking  ratios.  However,  this  is  not  expected  to  affect  the 
initial  rate  of  stacking  or  the  empirically  based  expectation  that  the  dye  concentration  should 
level  off  after  ~100  s.  This  simulation  shows  no  sign  of  such  a  steady  state.  To  investigate 
whether  the  steady  state  shows  up  at  higher  values  of  the  field  strength,  a  case  with  a  ten  times 
higher  field  strength  was  simulated,  as  described  in  the  next  section. 


Sample  Stacking  in  ID  at  40  kV/cm 

The  case  described  in  this  section  is  identical  to  the  one  in  the  previous  section,  except  for  the 
fact  that  the  channel  was  only  5  mm  long,  and  an  average  field  strength  of  2  MV/m  was 
described  (resulting  in  a  40kV/cm  effective  field  strength  over  the  2.5  mm  of  the  domain  that  has 
the  low  conductivity  buffer  (TE  zone)).  Using  a  time  step  of  1  ps,  this  case  was  simulated  for  a 
total  of  9  seconds. 
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Figure  4.13  shows  the  electrical  field  strength  for  this  case.  The  field  strength  starts  off  at  a  little 
over  40kV/cm  in  the  TE  zone,  and  increases  to  almost  60kV/cm  as  the  resistance  of  the  channel 
decreases  over  time.  With  this  higher  field  strength,  the  dye  concentration  increases  much  faster 
than  in  the  previous  case,  as  seen  in  Figure  4.14.  The  dye  begins  to  stack  immediately,  and 
continues  to  do  so  in  a  monotonic  fashion.  A  stacking  ratio  of  about  30,000  after  5  seconds,  and 
about  70,000  after  9  seconds  were  observed.  Note  that,  after  an  initial  strongly  nonlinear  stage, 
stacking  occurs  at  a  slightly  sublinear  pace.  Also,  after  9  sec,  no  steady-state  is  in  sight,  similar 
to  what  was  observed  in  the  previous  simulations. 
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Figure  4.13:  Evolution  of  the  field  strength. 

Figure  4.14:  Evolution  of  the  dye  eoneentration. 

Sample  Stacking  in  2D 

The  group  also  conducted  ITP  sample  stacking  studies  in  2D  under  limited  operating  conditions, 
corresponding  to  low  applied  voltage  difference  and  low  LE  molarity.  The  results  do  show 
stacking  of  the  dye  sample.  For  example,  with  applied  110  V/cm  and  a  0.01  M  NaCl  solution,  a 
stacking  ratio  of  30  after  2  sec,  and  200  after  20  sec  are  found.  However,  the  numerical 
construction,  which  is  robust  in  ID,  encountered  difficulties  as  stacking  progressed  in  the  2D 
context,  particularly  with  high  applied  voltage  difference  and  salt  concentration.  The  2D 
construction  is  ultimately  necessary  to  capture  dispersion  effects,  and  therefore  is  ultimately 
necessary  to  pursue.  This  is  left  as  a  matter  for  future  work. 
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4.2.6  Conclusions  and  Future  Work 


The  group  has  developed  and  implemented  a  detailed  eomputational  model  for  fundamental 
studies  of  isotachophoretie  sample  staeking  in  a  microchannel.  The  construction  accounts  for 
fluid  flow,  electrochemistry,  and  HEPES  buffer  physics,  including  non-electroneutral  bulk  flow 
conditions.  This  construction  was  used  to  conduct  parametric  studies  with  the  flow  geometry 
and  mixture  conditions  relevant  in  the  Stanford  experiments. 

Sample  stacking  behavior  using  ID  studies  were  observed  that  is,  in  many  respects,  in  good 
agreement  with  the  experimental  results.  The  detailed  structure  of  the  field  quantities  within  and 
outside  the  stacking  zone  is  generally  in  agreement  with  expected  behavior.  On  the  other  hand, 
this  agreement  is  mainly  qualitative.  Moreover,  even  qualitative  agreement  is  limited  to  the  early 
phase  of  stacking.  One-dimensional  results  clearly  do  not  replicate  the  observed  plateau  in 
stacking  ratio  seen  in  the  experiments.  Moreover,  a  rate  of  increase  of  the  stacking  ratio  that  is 
quantitatively  less  than  that  observed  experimentally  were  computed. 

Given  the  extent  of  detail  in  the  physical  models  implemented,  it  can  be  concluded  that  the  above 
disagreements  are  mainly  due  to  the  ID  nature  of  these  computational  results.  Long-term  2D 
stacking  studies  at  the  relevant  operating  conditions  remain  a  challenge,  and  are  a  subject  for 
future  work. 


5.0  SUMMARY  OF  PROJECT  IMPACT 

This  project  had  a  very  significant  impact  on  the  microfluidics,  lab-on-a-chip  and  biosensor 
communities  as  a  whole.  In  all,  the  project  resulted  in  over  168  journal  and  conference  papers, 
three  inventions  and  patents  (or  patent  applications),  and  trained  10  PhD’s,  three  MS  graduates, 
and  helped  trained  six  postdocs.  The  work  funded  as  part  of  this  project  also  led  to  six  awards 
from  conferences,  national  organizations,  and  the  US  government.  In  addition,  this  project 
funding  and,  more  importantly,  the  associated  technical  accomplishments  from  this  project 
contributed  to  the  successful  promotion  to  Associate  Professor  with  tenure  for  two  of  the  co-PIs 
involved  in  the  effort  (Myszka  and  Santiago).  Three  of  these  PhDs  and  postdocs  are  now  faculty 
members  at  major  universities  and  building  a  career  on  microfluidics  research.  Nine  others  are 
engineers,  researchers,  and/or  managers  concentrating  on  microfluidics  research  and 
development  at  either  companies  or  national  laboratories. 

The  models  produced  also  resulted  in  significant  technology  transfer  in  the  way  of 
commercialized  multiphysics  code  for  the  microfluidics.  Our  success  in  developing  the 
fundamentals  to  extract  accurate  information  from  protein  arrays  led  one  of  the  co-PIs  (Myszka) 
to  help  establish  a  micro  fluidic-based  company  in  Utah  called  Wasatch  Microfluidics. 

This  section  summarizes  each  of  these  accomplishments  and  ends  with  a  complete  listing  of  the 
50  journal  papers  which  resulted  from  this  project. 
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5.1  Training  of  Future  Leaders  in  Microfluidics 

The  following  students  obtained  their  MS  (as  their  terminal  degree)  largely  as  a  result  of  the 
work  they  eonducted  as  part  of  this  project: 

MS  Graduates  from  this  Project: 

•  Brian  Kidd,  M.S.  Meeh.  Eng.,  Stanford,  April  2001  (now  a  Researcher,  Stanford  Medical 
School) 

•  Bruce  P.  Mosier,  Eng.  Deg,  Meeh.  Eng.,  May  2001  (now  a  Research  Scientist  at  Sandia 
National  Labs,  working  on  Microfluidics) 

PhD  Graduates  from  this  Project: 

The  students  below  obtained  their  PhD  (as  their  terminal  degree)  largely  as  a  result  of  the  work 
they  conducted  as  part  of  this  project.  Also  listed  here  is  the  PhD  thesis  title  for  each  student. 

•  Chuan-Hua  Chen,  Ph.D.,  Meeh.  Eng.,  Stanford,  Jan.  2004,  Thesis  Title:  Microscale 
Electronkinetic  Transport  and  Stability,  Jan.  2004  (now  a  Research  Scientist  at  Rockwell 
Scientific  working  on  microfluidics) 

•  Amy  Herr,  Ph.D..  Meeh.  Eng,  Sept.  2002,  Title:  Isoelectric  Focusing  for  Multi-Dimensional 
Separations  in  Microfluidic  Devices,  Aug.  2002(currently  a  Research  Scientist  in  microfluidics 
at  Sandia  National  Labs;  to  start  a  faculty  position  at  University  of  California  at  Berkeley 
working  in  Microfluidics) 

•  Josh  I.  Molho, ,  Ph.D.  Meeh.  Eng.,  Stanford,  Dec.  2001,  Thesis  Title:  Electrokinetic 
Dispersion  in  Microfluidic  Separation  Systems(R&D  Engineer  and  Manager  at  Caliper 
Technologies;  a  microfluidics  company) 

•  Shankar  Devasenathipathy,  Thesis  Title:  Particle  Imaging  Diagnostics  and  Stacking 
Dynamics  in  Microfluidic  Systems,  Dec.  2003  (Researcher  Intel  Corporation;  working  on 
microfluidics  for  electronics  cooling) 

•  Michael  Oddy,  Ph.D.  Meeh.  Eng.,  Stanford,  December  2004,  Thesis  Title:  Electrokinetic 
Transport  Phenomena:  Mobility  Measurement  and  Electrokinetic  Instability,  Feb.  2005  (Finance 
Quantitative  Analyst,  Barclays  Capital) 

•  Rajiv  Bharadwaj,  Ph.D.  Chem.  Eng.,  Stanford,  June  2005,  Thesis  Title:  Microscale 
Electrokinetic  Sample  Stacking  (currently  R&D  Engineer,  Caliper  Technologies;  a 
microfluidics  company) 
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•  Shuhuai  Yao,  Ph.D.  Mech.  Eng.,  Stanford,  Sept.  2005,  Thesis  Title:  Theory,  Design,  and 
Demonstration  of  Electroosmotic  Pump  Technologies,  Aug.  2005(Postodoctoral  researcher  at 
Sandia  National  Labs  working  on  microfluidic  devices  for  protein  folding  studies) 

•  Dahzi  Wang,  PhD  2003,  University  of  California  at  Santa  Barbara,  Thesis,“  Micro  PIV 
Measurements  in  AC  Electrothermal  Flow”. 

Postdoctoral  Researchers  from  this  Project: 

•  Jim  Mikkelson,  Ph.D.  Chemistry,  Brown  Univ.,  left  Stanford  June  2002.  (Scientist  at 
Caliper  Technologies  working  on  microfluidics.) 

•  Peng  Zhou,  Ph.D.  Mech.  Eng.,  Stanford,  left  Stanford  September  2003.  (R&D  Engineer, 
Cooligy,  Inc.;  a  microfluidics  company  co-founded  by  Santiago.) 

•  Hao  Lin,  Ph.D.,  University  of  California  at  Berkeley,  left  Stanford  June  2005.  (Currently  a 
tenure-track  professor  at  Rutgers  University  working  on  microfluidics.) 

•  Jonathan  D.  Posner,  Ph.D.,  Univ.  of  California  at  Irvine,  left  Stanford  September  2005. 
(Currently  a  Professor,  Arizona  State  University  working  on  microfluidics.) 

•  Yasmina  Abdiche,  Ph.D.,  Biochemistry  Oxford  Univ.,  left  Utah  March  2005.  (Currently  a 
Senior  Research  Scientist  at  Pfizer  Inc.,  applying  expertise  in  biosensor  technology  to  drug 
discovery.) 

•  Cheryl  Baird,  Ph.D.,  Biochemisty  Univ.  of  Utah,  left  Utah  in  December  2004  to  work  on  array- 
based  biosensor  technology  at  HTS  Biosystems  (Boston,  MA).  Currently  works  at  Pacific 
Northwest  National  Laboratory  on  protein  microarrays  and  novel  detection  platforms. 

As  mentioned  above,  co-PIs  Myszka  and  Santiago  also  leveraged  the  technical  accomplishments, 
publications,  and  funding  from  this  project  to  earn  promotion  to  Associate  Professor  with  tenure 
in  their  respective  institutions.  Associate  Professor  and  co-PI  Meinhart  received  tenure  a  few 
months  after  notice  of  this  award  was  given. 

5.2  Technology  Transfer 

During  the  course  of  this  project,  several  simulation  models  have  been  developed  and 
implemented  in  CFD-ACE+.  Using  these  models,  several  physical  phenomena  have  been 
studied  in  detailed.  These  studies  directly  benefited  the  project  members,  as  well  as  academic 
and  industry  community  at  large.  These  accomplishments  are  discussed  in  detail  next: 


Models  Developed  (an  “M”  prefix  in  the  list  below  indicates  a  model  developed  as  part  of  this 
project): 
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Ml.  Particle  dielectrophoreis  models  that  include  eonventional  DEP,  traveling  wave  DEP, 
and  electrorotation  [29] 

M2.  AC  electrohydrodynamic  phenomenon  of  eleetrothermal  flows  [16] 

M3.  DC  electroosmosis  in  systems  with  thiek  double  layer  [P23] 

M4.  Slip  boundary  models  for  flow  in  hydrophobic  microchannel  systems  [33,  34] 

M5.  Models  for  EOF  mobility  variations  as  a  funetion  of  pH  for  fused  siliea  substrates[35] 

M6.  Unit-eell  based  models  for  highly  eoneentrated  buffer  solutions.  Models  prediet 
following  physical  properties[P23]: 

d.  Electrophoretic  mobility 

e.  Diffusion  coefficient 
f  Eleetrical  conduetivity 

M7.  Integration  of  biochemistry  models  with  partiele  DEP  models  [P45] 

M8.  Joule  heating[36] 

M9.  Super  LU  solver  for  highly  non-net  neutral  system[37] 

MIO.  Generalized  surface  ehemistry  models  (beads  and  surfaces)[38] 

Mil.  Biochemistry  database  integrated  with  electrochemistry 

M12.  Hydrogel  models  for  surface  biochemistry [30] 

Ml  3.  Least  square-based  engine  for  extraetion  of  kinetic  coefficients [3 8] 

Ml 4.  Rapid  ANN  based  kinetie  coeffieient  extraetion  algorithm  (this  has  not  been 

published  yet) 

Ml 5.  Python  based  scripts  to  automate  geometric  generation  for  proteomie  ehips 

(these  scripts  are  available  from  ESI-CFD,  NA,  Huntsville,  AL) 

The  above  models  been  validated  with  following  analytieal  solutions  (“S”  indicates  an  analytieal 
solution): 

5 1 .  Joule  heating  in  parallel  plate 

52.  Partiele  DEP  in  periodic  electrode  array 

53.  Partiele  DEP  in  wedge  shaped  (conical)  electrode  system 

54.  Eleetrothermal  flow  in  periodie  eleetrode  array 

55.  Ion  transport  near  charged  surfaee 

56.  EOF  flow  in  arbitrary  zeta  potential  and  double  layer  thiekness 

57.  All  surfaee  ehemistry  models 
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58.  After  validation,  the  models  have  been  used  extensively  to  study  the  following 
phenomena: 

59.  Sample  dispersion  under  non-uniform  zeta  potential 

510.  Eleetrothermal  flow  in  various  eleetrode  eonfigurations  and  elueidate  relative 
significanee  of  eleetrothermal  and  DEP  forces 

511.  Field  Amplified  Sample  Stacking 

51 2.  Effect  of  slip  at  nano-scale  on  overall  pressure  drop  in  microfluidic  systems 

513.  Manipulation  of  particles  in  AC  electric  fields  using  traveling  wave  DEP 

514.  Influence  of  double  layer  thickness,  zeta  potential  and  surface  charges  on  EOF  at  nano¬ 
scale 

SI  5.  Transport  of  analytes  in  highly  concentrated  buffer 

51 6.  Bead  based  immunoassays 

51 7.  Mass  transport  limitations  on  sensor  performance 

SI  8.  Effect  of  hydrogel  layer  on  proteomic  chip  performance 

5.2.1  Commercialization  of  Developed  Technologies 


Based  on  the  knowledge  and  understanding  of  above  physical  phenomena,  the  team  has  designed 
and  optimized  various  microfluidic  devices  such  as: 

D 1 .  Traveling  wave  DEP  based  Flow  Field  Fractionation  device 

D2.  Sample  Stacking  device 

D3.  Tunable  Laser  Cavity  (TLCC)  sensor 

D4.  Nanofluidic  systems 

D5.  HTS  proteomic  chip 

D6.  BIACORE  biosensors 

D7.  Notional  immunoassay  based  sensors 

Design  guidelines  have  been  generated  for  engineers  and  scientists  who  work  in  this  field. 
Besides,  the  simulation  models  have  been  implemented  in  CFD-ACE+  Version  2002,  2003  and 
2004,  and  have  been  successfully  marketed. 

Commercial  Release  of  Developed  Software  Capabilities:  Both  the  electrokinetic  models  and 
the  biochemical  assay  models  have  also  been  fully  integrated  into  the  commercial  CFD-ACE+ 
code,  along  with  the  development  of  associated  GUI,  user  manuals  and  tutorial  test  cases. 

Master  Classes  for  Advanced  CFD-ACE+  Users:  Sample  test  cases  were  also  included  based  on 
these  capabilities  in  the  Master  Classes  (for  advanced  CFD-ACE+  users)  for  the  Biochemistry 
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and  the  electrokinetic  modules.  These  classes  were  taught  at  the  Annual  CFD-ACE+  Users 
Conference. 

Multi-physics  Courses:  This  material  (test  cases,  tutorial  examples,  etc.)  is  also  be  included  in 
the  following  two  Biotechnology  industry-specific  multiphysics  courses  being  developed  at 
CFDRC,  which  were  offered  starting  November,  2001  in  various  locations  around  the  country. 
The  two  courses  offered  dealt  with:  (a)  design  of  biochemical  assays  in  microfluidic  systems, 
and  (b)  electrokinetic  and  electrochemistry  applications  in  microfluidic  systems.  The  objective 
of  these  courses  was  to  facilitate  rapid  learning  and  confidence  building  in  existing  and  new 
users  of  the  current  software  by  discussing  theoretical  background  (adequate  physical  models, 
data  sources,  and  boundary  conditions),  optimal  computational  strategies,  model  set-up  protocols 
of  CFD-ACE+,  and  expected  results  (accuracy,  speed,  uncertainty). 


5.2.2  Patents 

The  following  is  a  list  of  patents  and  patent  applications  resulting  from  this  project.  The  “I” 
prefix  in  the  numbers  below  indicate  an  invention  resulting  in  a  patent  application: 

11.  Santiago,  J.G.,  M.  Oddy,  and  J.C.  Mikkelsen,  2001.  “Electrokinetic  Instability 
Micromixer,”  Patent  Application  No.  20020125134.  September  12,  2002. 

12.  Santiago,  J.G.,  and  B.  Jung,  “A  Novel  Sample  Stacking  Capillary  Electrophoresis 
Device,”  Applied  for  provisional  patent,  2003. 

13.  Santiago  J.G.,  B.  Jung,  and  R.  Bharadwaj,  "Microfluidic  Sample  Stacking  Method  Using 
Single-Interface  Isotachophoresis",  Applied  for  provisional  patent,  2006. 

5.3  Awards  and  Honors 

The  following  awards  and  honors  resulted  at  least  in  part  from  the  research  conducted  and 
inventions  produced  as  part  of  this  project: 

Santiago  elected  Co-Chair  of  Gordon  Conference  on  the  Physics  and  Chemistry  of  Micro  fluidics, 
Oxford,  England  (planning  for  conference  in  August,  2007),  2005  (this  honor  is  a  general 
recognition  of  the  success  of  the  current  microfluidics  work) 

Santiago  group:  Best  Poster  in  session  award,  Gordon  Conference  on  the  Physics  and  Chemistry 
of  Microfluidics,  Oxford  England,  2005  (on  electrokinetic  instability) 

Santiago:  Presidential  Early  Career  Award  for  Scientist  and  Engineers  (PECASE),  2004  (This  is 
the  highest  award  by  the  US  Government  to  scientist  starting  their  career;  it  was  awarded 
for  the  Santiago  group’s  work  on  electrokinetic  instabilities  and  sample  stacking) 

Santiago  and  Huber,  Best  Poster  award  Annual  Meeting  of  the  American  Institute  of  Chemical 
Engineering  and  American  Electrophoresis  Society,  San  Francisco,  California,  2003  (on 
sample  stacking  work) 
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Santiago’s  group  received  a  Best  Poster  award  Gordon  Conference  on  the  Physics  and  Chemistry 
of  Microfluidics,  Big  Sky  Montana  2003  (on  sample  stacking  work) 

Santiago  and  M  Oddy  won  a  National  Inventors  Hall  of  Fame  Collegiate  Inventors  Award,  2001 
(on  the  design  of  an  electrokinetic  instability  micromixer) 

5.4  Complete  List  of  Publications  Resulting  from  this  Project 

The  “P”  prefix  in  the  numbered  list  below  indicates  a  publication  (journal  paper  or  conference 
paper)  resulting  from  this  project: 

Stanford  Publications  and  Presentations 

PI.  Jung,  B.,  Zhu,  Y.,  and  J.G.  Santiago,  "Detection  of  100  attomolar  Fluorophores  Using  a 
High  Sensitivity  On-Chip  CE  System  and  Transient  Isotachophoresis,"  submitted  to 
Analytical  Chemistry,  2006. 

P2.  Jung,  B.,  Bharadwaj,  R.,  Santiago,  J.G.,  "On-Chip  Million-Fold  Sample  Stacking  Using 
Single-Interface  Isotachophoresis,"  Vol.  78,  No.  7,  Analytical  Chemistry,  pp.  2319- 
2327,  2006. 

P3.  Bharadwaj,  R.  and  J.G.  Santiago,  "Dynamics  of  Field  Amplified  Sample  Stacking," 
Journal  of  Fluid  Mechanics,  Vol.  543,  57-92,  2005. 

P4.  Lin,  H.,  Storey,  B.,  M.  Oddy,  Chen,  C.-H.,  and  J.G.  Santiago,  “Instability  of 

Electrokinetic  MicroChannel  Flows  with  Conductivity  Gradients,”  Physics  of  Fluids, 
Vol.  16,  No.  6,  2004 

P5.  Alexis- Alexandre,  G.,  B.  Mohammadi,  J.G.  Santiago,  and  R.  Bharadwaj,  “Microfluidic 
Flow  Simulations:  Stacking  One-Dimensional  Study,”  Houille  Blance-Revue 
Internationale  De  Leau,  No.  5,  pp.  18-23,  2003. 

P6.  Jung,  B.,  R.  Bharadwaj,  and  J.G.  Santiago,  “Thousand-Fold  Signal  Increase  using  Field 
Amplified  Sample  Stacking  for  On-Chip  Electrophoresis,”  Electrophoresis,  Vol.  24, 
No.  19-20,  pp.  3476-3483,  2003 

P7.  Herr,  A.E.,  J.I.  Molho,  K.A.  Drouvalakis,  J.C.  Mikkelsen,  P.J.  Utz,  J.G.  Santiago,  and 
T.W.  Kenny,  “On-Chip  Coupling  of  Isoelectric  Focusing  and  Free  Solution 
Electrophoresis  for  Multi-Dimensional  Separations,”  Nna/yt/ca/  Chemistry,  Vol.  75, 
No.  5,pp.  1180-1187,  2003. 

P8.  Devasenathipathy,  S.,  J.G.  Santiago,  S.T.  Wereley,  and  C.D.  Meinhart,  "Particle  Tracking 
Techniques  for  Microfabricated  Fluidic  Systems,"  Experiments  in  Fluids,  Vol.  34, 

No.  4,pp.  504-513,2003 

P9.  Chen,  C.-H.,  Lin,  H.,  Lele,  S.K.  &  Santiago,  J.G.,  "Convective  and  absolute 

electrokinetic  instability  with  conductivity  gradients".  Journal  of  Fluid  Mechanics, 
524,  pp.  263  -  303,  2005. 

PIO.  Bharadwaj,  R.,  J.G.  Santiago,  and  B.  Mohammadi,  "Design  and  Optimization  of  On- 
Chip  Capillary  Electrophoresis,"  Electrophoresis,  Vol.  23,  pp.  2729-2744,  2002 

PIT  Mosier,  B.P.,  J.I.  Molho,  and  J.G.  Santiago,  "Bleached-Fluorescence  Imaging  of 
Microflows,"  Experiments  in  Fluids,  Vol.  33,  No.  4,  pp.  545-554,  2002 
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P12.  Devasenathipathy,  S.,  J.G.  Santiago,  and  K.  Takehara,  "  Particle  Tracking  Techniques 
for  Electrokinetic  MicroChannel  Flows,"  Analytical  Chemistry,  Vol.  74,  No.  15,  pp. 
3704-3713,  2002 

P13.  Chen,  C.H.,  and  J.G.  Santiago,  “Electrokinetic  Instability  in  High  Concentration 
Gradient  Microflows”,  ASME  Intemation  Mechanical  Engineering  Congress  & 
Exposition,  MEMS  and  Nanotechnology  Symposium,  November  17-22,  2002,  New 
Orleans. 

P14.  Devasenathipathy,  S.,  R.  Bharadwaj,  and  J.G.  Santiago,  “Investigation  of  field 

amplified  sample  stacking  with  particle  image  velocimetry,”  ASME  International 
Mechanical  Engineering  Congress  &  Exposition,  MEMS  and  Nanotechnology 
Symposium,  November  17-22,  2002,  New  Orleans. 

P15.  Bharadwaj,  R.,  J.G.  Santiago,  and  B.  Mohammadi,  "Design  and  Optimization  of  On- 
Chip  Capillary  Electrophoresis,"  Electrophoresis,  Vol.  23,  pp.  2729-2744,  2002. 

PI 6.  Devasenathipathy,  S.  and  J.G.  Santiago,  "Electrokinetic  Flow  Diagnostics,"  in  press, 
Microdiagnostics,  New  York,  Springer  Verlag,  2002. 

PI 7.  Devasenathipathy,  S.,  J.G.  Santiago,  and  K.  Takehara, "  Particle  Tracking  Techniques 
for  Electrokinetic  MicroChannel  Flows,"  Vol.  74,  No.  15,  pp.  3704-3713,  Analytical 
Chemistry,  2002. 

P18.  Mosier,  B.M.,  J.I.  Molho,  and  J.G.  Santiago,  "Bleached-Fluorescence  Imaging  of 
Microflows,"  in  press,  Experiments  in  Fluids,  2002. 

P19.  Mohammadi,  B.,  J.I.  Molho,  J.G.  Santiago,  “Incomplete  Sensitivities  in  the  Design  of 
Minimal  Dispersion  Fluidic  Channels,”  in  press,  Computational  Methods  in  Applied 
Mechanics  and  Engineering,  2002. 

P20.  Devasenathipathy,  S.,  J.G.  Santiago,  S.T.  Wereley,  and  C.D.  Meinhart,  "Particle 
Tracking  Techniques  for  Microfabricated  Fluidic  Systems,"  submitted  to  the 
Experiments  in  Fluids,  2001. 

P21.  Oddy,  M.H.,  J.G.  Santiago,  and  J.C.  Mikkelsen,  Jr.,  "Electrokinetic  Instability 

Micromixing,"  Vol.  73,  Electrokinetic  Instability  Micromixing,  pp.  5822-5832,  2001. 

P22.  Mohammadi,  B.,  J.G.  Santiago,  “Simulation  and  design  of  extraction  and  separation 
fluidic  devices",  Mathematical  Modelling  and  Numerical  Analysis,  Vol.  35,  No.  3, 
pp.  513-523,2001. 

P23.  Santiago,  J.G.,  "Electroosmotic  Flows  in  Microchannels  with  Finite  Inertial  and 
Pressure  Forces,"  Analytical  Chemistry,  Vol.  73,  No.  10,  pp.  2353-2365,  2001. 

P24.  Bharadwaj,  R.,  B.  Mohammadi,  and  J.G.  Santiago,  "Incomplete  Sensitivities  in 
Design  and  Control  of  Fluidic  Channels,"  Center  for  Turbulence  Research,  Stanford 
University,  Annual  Research  Briefs,  2001. 

P25.  Bharadwaj,  R.  and  J.G.  Santiago,  "Dynamics  of  Field  Amplified  Sample  Stacking," 
Proceedings  of  the  International  Mechanical  Engineering  Congress  and  Exposition, 
Sixth  Micro-Fluidic  Symposium,  New  York,  New  York,  2001. 
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P26.  Oddy,  M.H.,  J.G.  Santiago,  and  J.C.  Mikkelsen,  Jr.,  "An  Electrokinetic  Instability 
Micromixer,"  Proceedings  of  the  International  Mechanical  Engineering  Congress  and 
Exposition,  Sixth  Micro-Fluidic  Symposium,  New  York,  New  York, 
IMECE2001/MEMS-23882,  2001. 

P27.  Herr,  A.E.,  J.I.  Molho,  R.  Bharadwaj,  J.G.  Santiago,  and  T.W.,  Kenny,  "Miniaturized 
Isoelectric  Focusing  as  a  Component  of  a  Multi-dimensional  Separation  System," 
Fifth  International  Symposium  on  Micro  Total  Analysis  Systems  (pTAS)  Monterey, 
California,  pp.  51-53,  2001. 

P28.  Bharadwaj,  R.  and  J.G.  Santiago,  "Optimization  of  Field  Amplified  Sample  Stacking 
on  a  Microchip,"  Fifth  International  Symposium  on  Micro  Total  Analysis  Systems 
(pTAS)  Monterey,  California,  pp.  613-614,  2001. 

P29.  Oddy,  M.H.,  J.G.  Santiago,  and  J.C.  Mikkelsen,  Jr.,  "Electrokinetic  Instability 
Micromixing,"  Fifth  International  Symposium  on  Micro  Total  Analysis  Systems 
(pTAS)  Monterey,  California,  2001. 

P30.  Chen,  C.-H.,  J.C.  Mikkelsen,  Jr.,  J.G.  Santiago,  "Electrophoretic  Band  Crossing  for 
Measurements  of  Biomolecular  Binding  Kinetics,"  Technical  Digest  of  International 
Forum  on  Biochip  Technologies, ,  2000.Beijing,  China,  pp.  441-442,  Oct  11-14  2000. 

P31.  Herr,  A.E.,  J.C.  Mikkelsen,  Jr.,  J.G.  Santiago,  T.W.,  Kenny,  "Electroosmotic  Flow 
Suppression  and  its  Implications  for  A  Miniaturized  Full-Field  Detection  Approach  to 
Capillary  Isoelectric  Focusing  (cIEF),"  Proceedings  of  ASME  2001  Winter  Annual 
Meeting,  Fifth  Micro-Fluidic  Symposium,  International  Mechanical  Engineering 
Congress  and  Exposition,  Orlando,  Florida,  MEMS  Vol.  2,  pp.  513-518,  2000. 

P32.  Mosier,  B.P.,  J.C.  Mikkelsen,  Jr.,  J.G.  Santiago,  "  A  Novel  Bleached  Fluorescence 
Imaging  Technique  for  Microfluidics,"  Fifth  Micro-Fluidic  Symposium,  International 
Mechanical  Engineering  Congress  and  Exposition,  Orlando,  Florida,  2000. 

P33.  Devasenathipathy,  S.,  J.I.  Molho,  J.C.  Mikkelsen  Jr.,  J.G.  Santiago,  "Electroosmotic 
Flow  Field  Measurements  with  Particle  Image  Velocitry,"  Fifth  Micro-Fluidic 
Symposium,  International  Mechanical  Engineering  Congress  and  Exposition, 
Orlando,  Florida,  2000. 

P34.  Herr,  A.E.,  J.I.  Molho,  J.G.  Santiago,  T.W.  Kenny, D.  A.  Borkholder,  G.J.  Kintz,  P. 
Belgrader,  M.  A.  Northrup,  "A  Miniaturized  Full-Field  Array  Detection  Approach  to 
Protein  Separations  Using  Capillary  Isoelectric  Focusing  (cIEF),"  Proceedings  of  the 
Micro  Total  Analysis  Systems  Symposium,  Enschede,  The  Netherlands,  pp.  367-370, 
2000. 
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CFDRC  Corporation  Publications  and  Presentations 


P35.  Przekwas,  A.  J.  and  V.  B.  Makhijani,  “Mixed-Dimensionality,  Multi-Physics 
Simulation  Tools  for  Design  Analysis  of  Microfluidic  Devices  and  Integrated 
Systems”,  In:  Technical  Proc.  4th  International  Conf  On  Modeling  and  Simulation  of 
Microsystems  (MSM  2001),  Hilton  Head,  SC,  Laudon,  M  and  Romanowicz,  B, 
(Eds.),  Computational  Publications,  Cambridge,  MA,  pp.  198-201,  2001. 

P36.  V.  B.  Makhijani  and  S.  Sundaram,  “High-fidelity  Computational  Modeling  and 
Analysis  of  Sensors  for  Biowarfare  Agent  Detection”,  presented  at  Lockheed  Martin 
NESS,  Manassas,  VA,  May  2,  2001. 

P37.  A.  K.  Singhal  and  V.  B.  Makhijani,  “Growing  Role  of  Multi-Physics  Simulations  for 
Rapid  Development  of  Micro-Reaction  and  Bio-Analytical  Systems”,  In  Proc.  5* 
International  Conf  on  Microreaction  Technology  (IMRET5),  Strasbourg,  France, 
May  27-30,  2001. 

P38.  S.  Krishnamoorthy,  J.J.  Feng  &  V.  B.  Makhijani,  “Full-Scale  Computational  Study 
Of  Electrophoretic  Systems”,  In.  Proc  Electrochemical  Society  2001  Joint 
International  Meeting  San  Fransisco,  CA,  Sept.  4,  2001. 

P39.  S.  Krishnamoorthy  and  M.  Shariati,  “CFD-ACE-I-:  Multi-Physics  Simulation  Software 
for  Computer-Aided  Design  and  Optimization  of  Biochip  Systems”,  Presented  at 
Caliper  Technologies  Corporation,  Mountain  View,  CA,  Sept.  5,  2001. 

P40.  S.  Krishnamoorthy  and  M.  Shariati,  “CFD-ACE+:  Multi-Physics  Simulation  Software 
for  Computer-Aided  Design  and  Optimization  of  Biochip  Systems”,  Presented  at 
Lynx  Therapeutics,  Inc.,  Hayward,  CA,  Sept.  6,  2001. 

P41.  J.  Feng,  S.  Krishnamoorthy  and  V.  B.  Makhijani,  "Simulation  of  Dielectrophoresis  of 
Particles  Near  Microelectrodes",  In:  Proc.  Labautomation  2002,  Palm  Springs,  CA, 
Jan.  27  -29,  2002. 

P42.  Feng,  J.  J.,  Blosch,  E.,  Krishnamoorthy,  S.  and  Makhijani,  V.  B.,  “Computational 
Analysis  of  Electrokinetic  Transport  in  Lab-On-a-Chip  Systems”,  In:  Proc.  IBC’s  8*'^ 
Annual  International  Microtechnology  Event.  Chips  to  Hits  ’01.  San  Diego.  CA.  Oct. 
28* -Nov.  1st,  2001. 

P43.  Prabhakarpandian,  B.,  Jenkins,  J.,  Sundaram,  S.  and  Makhijani,  V.  B.,  “Multi-Physics 
Simulations  For  The  Design  of  Novel  and  Optimal  Biosensors  &  Biodiagnostics 
Systems”,  In:  Proc.  IBC’s  8*  Annual  International  Microtechnology  Event,  Chips  to 
Hits  ’01,  San  Diego,  CA,  Oct.  28* -Nov.  1st,  2001. 

P44.  Krishnamoorthy,  S.,  Feng,  J.  J.,  Chen,  Z.  J.  and  Makhijani,  V.  B.,  “Numerical  and 
Analytical  Studies  of  AC  Electric  Field  in  Dielectrophoretic  Electrode  Arrays”,  In: 
Proc.  5*  International  Conference  on  Modeling  and  Simulation  of  Microsystems 
(MSM  2002),  San  Juan,  Puerto  Rico,  April  22-25,  2002. 

P45.  Prabhakarpandian,  B.,  Shah,  K.  B.,  Sundaram,  S.  and  Makhijani,  V.  B.,  “Biochemical 
Binding  in  Microspheres-Based  Immunoassays”,  In:  Proc.  5*  International 
Conference  on  Modeling  and  Simulation  of  Microsystems  (MSM  2002),  San  Juan, 
Puerto  Rico,  April  22-25,  2002.V.B.Makhijani,  S. Sundaram,  B.P.Pandian,  and  S. 
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